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ABSTRACT 


Spacecraft  with  high  performance  attitude  control  systems  requirements  have 
traditionally  relied  on  imperfect  mechanical  gyroscopes  for  primary  attitude 
determination.  Gyro  bias  errors  are  corrected  with  a  Kalman  filter  algorithm  that  uses 
updates  from  precise  attitude  sensors  like  star  trackers.  Gyroscopes,  however,  have  a 
tendency  to  degrade  or  fail  on  orbit,  becoming  a  life  limiting  factor  for  many  satellites. 
When  errors  become  erratic,  pointing  accuracy  may  be  lost  during  short  star  gaps. 
Unpredictable  gyro  degradations  have  impacted  NASA  spacecraft  missions  such  as 
Sky  lab  and  Hubble  Space  Telescope  as  well  as  several  DoD  and  ESA  satellites.  An 
alternative  source  of  angular  rate  information  is  a  software  implemented  real  time 
dynamic  model.  Inputs  to  the  model  from  internal  sensors  and  known  spacecraft 
parameters  enable  the  tracking  of  total  system  angular  momentum  from  which  body  rates 
can  be  determined.  With  this  technique,  the  Kalman  filter  algorithm  provides  error 
corrections  to  the  dynamic  model.  The  accuracy  of  internal  sensors  and  input  parameters 
determine  the  effectiveness  of  this  angular  rate  estimation  technique.  This  thesis  presents 
the  background  for  understanding  and  implementation  of  this  technique  into  a 
representative  attitude  determination  system.  The  system  is  incorporated  into  an  attitude 
simulation  model  developed  in  SIMULINK  to  evaluate  the  effects  of  dynamic  modeling 
errors  and  sensor  inaccuracies.  Results  are  presented  that  indicate  that  real  time  dynamic 
modeling  is  an  effective  method  of  angular  rate  determination  for  maneuvering  multi¬ 
body  spacecraft  attitude  control  systems. 
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I. 


INTRODUCTION 


A.  ANGULAR  RATE  ESTIMATION 

The  attitude  control  systems  of  today’s  high  tech  satellites  require  accurate 
angular  rate  knowledge  for  attitude  propagation  and  control  loop  damping  in  order  to 
meet  pointing  and  tracking  requirements.  Gyroscopes  offer  varying  degrees  of  precision 
for  direct  measurement  of  angular  rates  and  have  been  the  primary  attitude  determination 
sensor  used  by  spacecraft  for  many  years.  However,  high  cost  and  low  reliability  has 
lead  users  and  designers  to  explore  other  options.  Using  dynamic  calculations, 
uncertainty  based  estimation  algorithms  or  a  combination  of  these  two  methods,  onboard 
processors  can  estimate  spacecraft  angular  rates  without  measuring  them  directly. 
Angular  rate  estimation  techniques  can  be  a  viable  alternative  for  back-up  or  even 
primary  attitude  determination  system  in  many  control  schemes. 

B.  OPTIMAL  ESTIMATION 

“An  optimal  estimator  is  a  computational  algorithm  that  processes  measurements 
to  deduce  a  minimum  error  estimate,  in  accordance  with  stated  optimality  criterion,  of  the 
state  of  a  system  by  utilizing:  knowledge  of  system  measurement  dynamics,  assumed 
statistics  of  system  noises  and  measurement  errors,  and  initial  condition  information.” 
[Ref.  1]  Estimation  techniques  provide  filtering,  smoothing  and  prediction  of  state 
variables  in  an  imperfect  model  based  on  imperfect  measurement  update  data.  The  most 
common  optimal  estimator  used  in  stochastic  systems  is  the  Kalman  filter.  The  dynamic 
model  error  and  measurement  error  are  assumed  as  zero  mean  Guassian  white  noise 
processes  with  known  covariance.  Estimators  are  commonly  used  even  in  systems  where 
all  state  variables  required  by  controller  can  be  measured. 
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Maybeck  identifies  three  basic  reasons  why  deterministic  methods  are  insufficient 
for  describing  real  systems.  First,  no  mathematical  model  is  perfect.  There  are  always 
effects  that  are  necessarily  neglected  to  make  the  model  practical  and  even  modeled 
effects  are  only  approximations  to  what  is  physically  occurring.  Second,  dynamic 
systems  are  driven  not  only  by  control  inputs  but  also  by  disturbances  which  we  can 
neither  control  nor  model  deterministically.  Lastly  sensors  do  not  provide  perfect  and 
complete  data.  [Ref.  2] 

An  important  advantage  of  using  an  optimal  estimator  such  as  the  Kalman  filter  is 
that  the  attitude  determination  output  does  not  affect  the  controller  design.  Therefore  the 
development  of  an  optimal  controller  can  be  accomplished  independently.  An  attitude 
control  algorithm  optimized  for  ideal  state  inputs  will  remain  optimized  with  state 
estimates  provided  by  an  optimal  estimator. 

C.  ANGULAR  RATE  ESTIMATION  FROM  VECTOR  MEASUREMENTS 

The  concept  of  using  attitude  sensor  data  to  produce  an  estimate  of  spacecraft 
angular  rate  without  gyroscopes  is  not  new.  Gyroless  attitude  and  angular  rate  estimation 
has  been  a  prime  concern  for  small  inexpensive  spacecraft  that  do  not  carry  gyroscopes 
but  still  require  rate  information  for  attitude  propagation  and  control.  Estimation 
techniques  also  provide  options  for  complex  spacecraft  that  require  back-up  control 
modes  in  the  event  of  gyro  failures.  The  problem  of  angular  rate  estimation  can  be 


2 


treated  separately  from  attitude  determination.  Several  reliable  algorithms  have  been 
developed  that  produce  angular  rate  using  on-board  processors. 

It  is  possible  to  extract  angular  rate  directly  from  time  derivatives  of  measured 
vectors  resolved  in  the  body  coordinate  frame  and  known  in  an  inertial  frame  from  a 
model  or  almanac  data.  This  technique,  however,  requires  at  least  two  non-parallel 
vector  measurements.  It  also  exhibits  time  lag  and  produces  very  noisy  data  since  the 
algorithms  depend  on  derivatives  of  already  noisy  measurements.  A  better  estimate  of 
angular  rate  can  be  achieved  by  treating  the  problem  as  stochastic.  For  this  method  an 
estimating  filter  is  applied  that  uses  dynamic  equations  of  motion  to  take  advantage  of 
past  data  in  a  recursive  sense.  This  method  has  the  additional  advantage  of  being  able  to 
update  rate  estimates  at  time  steps  when  only  one  vector  measurement  is  available.  The 
vector  derivative  is  treated  as  a  noisy  measurement  update  to  a  Kalman  rate  estimator. 

A  Kalman  filter  requires  a  linear  state -space  dynamic  model  but  the  dynamics  of  a 
spacecraft  in  the  rotating  body  coordinate  frame  are  coupled  and  nonlinear.  Attitude 
determination  schemes  deal  with  this  problem  in  different  ways.  Several  methods  have 
been  proposed  for  interlacing  two  or  three  linear  Kalman  filters  together  to  capture  the 
nonlinear  dynamics  [Ref.  3].  It  is  shown  that  the  coupled  equations  of  motion  can  be 
written  as  a  linear  combination  of  the  spacecraft  angular  rate  and  two  other  newly  defined 
vectors  whose  components  are  nonlinear  combinations  of  the  angular  rate  vector 
elements.  The  differential  equations  for  each  of  these  new  vectors  can  be  written  as 
linear  combinations  of  the  other  and  the  angular  rate.  Adding  white  noise  vectors  to  these 
new  equations  of  motion  produces  a  set  of  three  stochastic  linear  models  that  can  be  used 
in  separate  Kalman  filters.  The  filters  are  interlaced  with  their  estimated  outputs  treated 
as  deterministic  inputs  to  each  other. 

An  effective  method  of  dealing  with  non-linear  dynamics  is  presented  for  use  by 
the  Pseudo-linear  Kalman  filter  (PSELIKA)  and  the  state -dependent  algebraic  Riccati 
equation  filter  (SDARE)  [Ref.  4].  For  these  unique  rate  estimation  techniques  the 
equation  of  motion  is  transformed  to  express  the  nonlinear  terms  in  angular  rate  as  a 
product  of  a  matrix  whose  elements  depend  on  the  components  of  the  angular  rate  vector 
and  the  angular  rate  vector  itself.  The  pseudolinear  Kalman  filter  or  PSELIKA  is 
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implemented  as  an  ordinary  discrete  Kalman  filter  with  a  time  dependent  state  transition 
matrix  whose  angular  rate  dependent  elements  are  formed  from  the  current  estimate  of 
angular  rate.  The  SDARE  or  state -dependent  algebraic  Riccati  equation  filter  is  also 
implemented  as  a  discrete  Kalman  filter  derived  from  the  same  representation  of  the 
dynamic  equation.  The  Kalman  gain  used  in  this  filter  is  computed  from  the  solution  of 
an  algebraic  Riccati  equation  involving  the  angular  rate  dependent  matrix.  This 
eliminates  the  need  to  propagate  and  update  the  filter  state  covariance  which  is  normally 
used  in  the  gain  calculation. 

The  pseudolinear  filter  concept  can  be  extended  to  use  quaternion  measurements. 
In  the  attitude  determination  system  on  board  the  Rossi  X-Ray  Timing  Explorer  (RXTE) 
satellite  advanced  star  tracker  software  directly  provides  information  in  the  form  of 
quaternions  [Ref.  5].  The  same  propagation  algorithm  is  used  with  measurement  updates 
in  the  form  of  quaternions.  Since  the  device  yielding  the  attitude  measurement  also 
conducts  a  star  search  and  identification  process  a  time  delay  is  introduced.  Two 
algorithms  are  presented  for  overcoming  this  delay  problem. 

D.  FULL  ORDER  ESTIMATORS 

It  is  not  necessary  to  treat  the  problem  of  angular  rate  estimation  separately  from 
attitude  determination.  In  fact,  attitude  sensor  data  also  requires  filtering  to  smooth  out 
measurement  noise  and  produce  a  clean  attitude  estimates.  High  precision  control  laws 
generally  require  both  attitude  and  angular  rate  information.  These  components  make  up 
the  full  system  state  vector. 

Full  order  estimators  produce  optimal  state  estimates  in  the  case  where  imperfect 
measurement  data  related  to  all  of  the  modeled  states  is  continuously  available. 
Spacecraft  with  a  triad  of  rate  gyroscopes  and  attitude  sensors  have  the  required  state 
measurements  but  the  bandwidth  of  the  attitude  sensors  is  normally  too  low  for  the 
control  system.  Additionally,  there  may  be  periodic  gaps  in  sensor  coverage  and  certain 
sensors  do  not  provide  reliable  attitude  information  in  all  three  orthogonal  dimensions. 
Therefore  standard  full  order  state  estimators  are  generally  not  used. 
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In  most  practical  implementations  of  gyro-based  attitude  control  systems  rate 
information  is  used  in  the  propagation  stage  of  an  attitude  estimator  based  on  kinematic 
equations.  The  estimator  utilizes  discrete  vector  observations,  resolved  in  both  body- 
fixed  coordinates  and  a  reference  frame,  to  estimate  spacecraft  attitude  and  gyro  drift 
rates. 

In  an  ideal  attitude  determination  system  where  perfect  knowledge  of  the 
spacecraft  angular  rates  is  available,  the  accuracy  of  the  kinematic  model  for  determining 
attitude  is  limited  only  by  the  processor  time  step  and  the  initial  attitude  state.  Attitude 
sensors  can  be  used  to  periodically  reinitialize  the  kinematic  model  but  the  data  they 
provide  is  discontinuous  and  corrupted  by  environmental  effects  and  sensor  design 
imperfections. 

In  real  attitude  determination  systems,  the  source  of  angular  rate  information  may 
be  noisy  as  well  as  biased,  as  is  the  case  with  spacecraft  gyroscopes.  This  results  in 
additional  errors  in  the  spacecraft  attitude  determination  system.  For  control  systems 
where  attitude  and  angular  rate  information  is  critical,  optimal  estimation  techniques  are 
employed  to  combat  the  effects  of  model  and  sensor  inaccuracies.  As  gyros  degrade  the 
estimator  can  be  made  more  robust  to  plant  error  at  the  expense  of  attitude  determination 
accuracy.  At  some  point  the  gyro  outputs  may  become  too  erratic  to  meet  control  system 
requirements. 

E.  REDUCED  ORDER  ESTIMATORS 

Obtaining  satellite  angular  rate  estimates  without  the  use  of  rate  gyroscopes  or 
other  deterministic  rate  data  can  be  accomplished  with  a  reduced  order  estimator.  This 
class  of  filters  produces  estimates  of  all  modeled  state  variables  when  only  a  subset  is 
directly  related  to  the  measurement  data.  In  this  case,  attitude  sensor  information  is 
available  but  no  direct  rate  measurement  is  performed.  Since  gyro  data  is  not  available 
no  update  to  the  spacecraft  rate  is  available  between  attitude  measurements. 

The  reduced  estimator  Kalman  filter  is  formulated  from  the  state  space  dynamic 
equations  of  motion  linearized  about  the  current  estimate  of  state.  Again,  for  spacecraft 
attitude  control  the  state  vector  includes  both  the  attitude  and  angular  rate.  The 
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deterministic  inputs  to  the  model  are  the  known  externally  applied  moments  usually 
consisting  only  of  control  torques.  Unmodeled  disturbances  and  other  dynamic  effects 
must  be  accounted  for  by  a  robust  plant  noise  model.  The  estimates  of  attitude  and 
angular  rate  are  propagated  forward  by  the  dynamic  model  with  discrete  corrections  from 
the  attitude  sensors. 

The  performance  of  the  reduced  Kalman  estimator  can  be  improved  by  including 
higher  order  dynamic  effects  in  the  system  plant  model  and  including  disturbance  torques 
as  deterministic  inputs.  However  this  costs  extra  processing  power,  requires  additional 
sensors  and  cannot  account  for  all  of  the  unknowns.  Additional  dynamic  complexities 
are  introduced  in  multi-body  satellites  that  have  time  varying  moments  of  inertia, 
changing  centers  of  mass  and  flexibility  modes.  Determining  the  magnitude  and  effects 
of  disturbances  and  modeling  simplifications  for  a  particular  spacecraft  is  an  important 
design  consideration  that  often  requires  rough  calculations,  simulation  and  engineering 
judgment.  In  general,  since  disturbances  are  of  low  amplitude  and  low  frequency 
compared  to  control  torques,  their  effects  can  be  accounted  for  as  plant  error  by  the  filter. 

One  approach  that  avoids  the  use  of  the  typically  uncertain  dynamic  model 
altogether  is  to  treat  the  spacecraft  as  a  noncooperative  target  whose  rotation  must  be 
tracked  by  a  stochastic  estimator  [Ref.  6].  This  concept  is  borrowed  from  tracking 
theory  where  it  has  been  widely  used  to  estimate  target  motion.  This  algorithm  adds 
angular  acceleration  to  the  state  vector  and  substitutes  attitude  states  with  the  integrated 
rate  parameters  to  formulate  a  nine  state  linear  Kalman  filter.  Since  this  set  of  variables 
serves  as  an  approximate  third-order  attitude  parameterization,  the  size  must  be 
controlled  by  sampling  interval.  Instead  of  using  dynamics,  time  propagation  of  the 
estimated  state  variables  is  performed  in  the  proposed  filter  by  modeling  the  spacecraft 
angular  acceleration  as  an  exponentially  autocorrelated  stochastic  process  and  using  a 
polynomial  kinematic  model. 

In  reduced  order  estimators,  time  steps  for  state  propagation  must  be  kept  small  to 
minimize  the  effects  of  dynamic  simplifications  and  linearization.  The  primary  limitation 
of  the  reduced  estimator  becomes  evident  during  prolonged  intervals  between  attitude 
sensor  measurement  updates.  Gaps  in  attitude  sensor  coverage  must  be  kept  short  enough 
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to  ensure  that  estimator  errors  do  not  grow  beyond  required  precision  of  the  attitude 
control  system. 

F.  STAR  SENSOR  BASED  ESTIMATORS 

Bandwidth  and  accuracy  limitations  of  attitude  sensors  have  precluded  the  use  of 
reduced  order  estimators  for  spacecraft  missions  with  high  performance  pointing  and 
tracking  requirements.  Aside  from  star  trackers,  most  sensors  do  not  provide  accurate 
and  stable  enough  attitude  data  for  high  performance  tracking  requirements. 
Additionally,  these  sensors  provide  discrete  measurements  unlike  gyroscopes  which 
produce  nearly  continuous  data  and  can  be  sampled  at  extremely  high  frequency  limited 
only  by  processing  capabilities.  However,  recent  advancements  in  sensor  technology 
suggest  that  star  trackers  can  provide  updates  accurately  and  frequently  enough  to  be 
considered  for  use  as  primary  attitude  sensors  for  a  wide  range  of  spacecraft  missions. 
Key  technological  improvements  include  a  wide  field  of  view,  high  sensitivity,  low  noise 
equivalent  angles  and  high  bandwidth  iteration  rate. 

It  has  been  proposed  that  advanced  high-bandwidth  star  sensors  can  be  used  as  the 
primary  sensors  for  on-orbit  attitude  rate  determination  in  place  of  angular  rate 
gyroscopes  [Ref.  7].  A  Kalman  filter  algorithm,  using  updates  from  star  tracker 
measurements,  is  presented  that  tracks  errors  in  the  attitude  and  angular  rate  for  a 
kinematic  based  state  space  model.  The  spacecraft  angular  rate  used  for  attitude 
propagation  is  separated  into  a  nominal  component  and  a  component  due  to  control  and 
disturbance  moments.  The  algorithm  assumes  that  the  bandwidths  of  the  control  system 
and  disturbance  effects  are  at  least  an  order  of  magnitude  smaller  than  the  measurement 
bandwidth.  This  allows  the  rate  of  change  in  spacecraft  angular  rate  due  to  unmodeled 
dynamic  effects  to  be  modeled  as  a  first-order  Gauss-Markov  process.  Nonlinear 
kinematics  are  used  in  the  attitude  propagation  phase  which  is  performed  external  to  the 
estimator.  The  errors  in  the  this  model  are  then  tracked  using  an  ordinary  linear  discrete 
Kalman  filter. 
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G.  ANALYTICAL  RATE  DETERMINATION 

Another  method  of  obtaining  estimates  of  spacecraft  angular  rates  is  through 
direct  calculation  from  the  dynamic  equations  of  motion  for  the  system.  This  rate 
calculation  is  performed  in  real  time  based  on  known,  derived  and  sensed  internal 
parameters  of  the  spacecraft.  The  software  implemented  dynamic  model  can  be  adjusted 
to  include  varying  levels  of  complexity  for  multi -body  spacecraft  that  operate  in  different 
configurations.  The  Aerospace  Corporation  has  patented  one  methodology  for  this 
process  called  the  “Pseudo  Gyro”  [Ref.  8]. 

Information  from  internal  sensors  that  detect  relative  orientations  and  angular 
rates  of  momentum  exchange  devices  and  appendages  are  used  to  continually  update  the 
parameters  used  in  the  dynamic  calculations.  These  parameters  determine  component 
contributions  to  the  system  angular  momentum  and  inertia  dyadic.  Control  torques  and 
modeled  disturbances  are  integrated  to  capture  external  dynamic  effects.  Using  the  total 
system  angular  momentum  and  inertia  dyadic  and  the  relative  momentum  of  internal 
components  and  appendages  from  the  dynamic  model  the  spacecraft  angular  rate  can  be 
calculated.  If  the  appendage  masses  are  small  or  their  relative  motions  are  slow  the 
spacecraft  inertia  matrix  is  slow  to  change.  The  moment  of  inertia  calculations  can  be 
performed  at  a  lower  frequency  to  save  processing  time.  Cross  product  effects  due  to  the 
rotating  coordinate  frame  must  also  be  accounted  for  in  the  dynamic  model  since  all 
measurements  are  taken  in  the  body  frame. 

The  accuracy  of  the  calculated  rate  is  dependent  on  the  quality  of  the  dynamic 
model  and  the  sensor  information  available.  Error  sources  include  imperfect  knowledge 
of  component  or  appendage  inertia  matricies  and  mass  centers  as  well  as  relative  angular 
position  and  rate  data  from  the  internal  sensors.  Often  relative  rates  are  not  measured 
directly  but  derived  from  position  encoders  which  adds  extra  noise  to  the  momentum 
along  the  axis  of  rotation  of  the  appendage.  Errors  are  also  introduced  through  external 
disturbance  effects  since  they  cannot  be  perfectly  modeled.  All  of  these  dynamic 
influences  effect  the  model  precision  in  varying  degrees.  Additionally,  due  to  the 
finite/discrete  processing  capabilities  of  the  computer  performing  the  calculations  there 
will  be  a  slight  drift  over  time  from  the  true  state  even  in  the  case  of  perfect  sensors  and 
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input  data.  This  drift  is  obviously  minimized  by  increasing  the  frequency  of  the  discrete 
model  calculations. 

H.  ATTITUDE  ESTIMATION  FROM  CALCULATED  RATE 

The  rate  estimate  produced  by  the  dynamic  calculations  can  be  fed  into  a  Kalman 
filter  in  place  measurement  data  from  a  gyroscope.  The  filter  receives  attitude  update 
information  from  attitude  sensors  to  produce  smoothed  attitude  reference  and  angular  rate 
error  estimates  to  enhance  satellite  attitude  determination  accuracy.  The  plant  error 
introduced  into  the  filter  is  a  combination  of  all  the  internal  sensor  and  modeling  errors 
that  are  used  as  inputs  to  the  dynamic  model.  All  known  sensor  biases  must  be 
incorporated  into  the  plant  model  since  the  Kalman  filter  assumes  errors  to  be  zero  mean 
Gaussian.  Depending  on  different  modes  of  operation,  different  plant  noise  models  may 
be  required  for  acceptable  filter  performance.  Normally,  filters  that  receive 
measurements  from  gyroscopes  are  designed  to  track  gyro  biases  which  remain  relatively 
constant  in  the  body  frame.  The  rate  error  from  the  dynamic  estimate,  however,  exhibits 
different  characteristics.  Since  the  dynamic  model  is  external  to  the  Kalman  filter  higher 
order  effects  can  be  more  easily  incorporated  into  the  dynamics  but  onboard  processing 
capability  may  still  limit  the  complexity. 

I.  PURPOSE 

The  objective  of  this  study  is  to  develop  and  evaluate  a  practical  attitude  and  rate 
estimation  scheme  for  a  multi-body  spacecraft  attitude  control  system  that  incorporates 
both  real  time  angular  rate  calculation  from  the  system  dynamic  model  (dynamic  gyro) 
and  a  Kalman  filter  estimator  with  attitude  sensor  updates  provided  by  star  trackers.  It  is 
hypothesized  that  acceptable  attitude  control  performance  can  be  realized  by 
maneuvering  multi-body  spacecraft  without  hardware  gyroscope  data  using  this 
methodology. 

The  evaluation  of  this  approach  is  performed  through  simulation  using  a  model 
developed  in  SIMULINK.  A  MATLAB  script  file  is  used  to  set  up  the  necessary 
initialization  parameters  and  system  specifications.  Simulation  results  are  presented 
graphically  in  MATLAB  plots.  The  performance  of  the  developed  gyroless  attitude 
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determination  system  is  compared  to  a  conventional  gyro-based  system  that  uses  the 
same  Kalman  filter  estimation  algorithm  and  attitude  updates.  Results  are  also  presented 
that  analyze  the  affects  of  several  major  error  sources  on  the  performance  of  this  dynamic 
gyro  based  attitude  determination  system. 
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II.  SIMULATION  OVERVIEW 


The  evaluation  and  analysis  of  attitude  determination  scheme  proposed  in  this 
thesis  is  done  through  modeling  and  simulation.  The  context  for  the  development  of  the 
attitude  determination  system  is  a  multi -body  attitude  control  system  for  a  maneuvering 
spacecraft.  An  attitude  simulation  model  for  the  spacecraft  that  includes  vehicle 
dynamics,  determination  and  control  subsystems  as  well  as  modeled  error  sources  is 
developed  using  MATLAB/SIMULINK.  An  overview  of  the  simulation  and  top  level 
subsystems  is  presented  in  this  chapter.  Subsequent  chapters  provide  the  descriptions  of 
the  subsystems  and  derivation  of  the  equations  on  which  the  model  is  based.  In  the 
following  chapters,  actual  SIMULINK  diagrams  are  presented  that  show  the 
implementation  of  associated  subsystems. 

A.  CONCEPTUAL  ATTITUDE  SIMULATION  MODEL 

A  conceptual  representation  of  the  simulation  developed  for  this  study  is 
illustrated  in  Figure  2. 
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Figure  2.  Conceptual  Attitude  Simulation  Block  Diagram 
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The  model  is  divided  into  subsystems  represented  in  color  shaded  blocks.  This 
breakdown  reduces  the  complexity  of  the  overall  model  into  manageable  segments  to  aid 
in  design  and  analysis.  Arrows  indicate  the  dynamic  coupling  and  flow  of  data  between 
subsystem  blocks.  The  top  level  SIMULINK  diagram  that  implements  the  concept  is 


shown  in  Figure  3. 


Figure  3.  Top  Level  SIMULINK  Attitude  Simulation  Model 


B.  SPACECRAFT  MOTION 

The  actual  spacecraft  attitude  motion  is  simulated  in  the  rotational  dynamics  and 
kinematics  subsystems  with  inputs  and  outputs  represented  by  solid  black  lines.  Relative 
motion  of  the  secondary  body  or  appendage  is  treated  in  a  separate  subsystem  with 
dynamic  effects  coupled  into  the  overall  spacecraft  motion  through  momentum  exchange 
directly  related  to  the  drive  motor  rate.  Dynamic  effects  of  reaction  wheel  control  are 
also  realized  through  momentum  exchange.  Magnetic  control  effects  are  input  to  the 
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spacecraft  dynamics  as  external  torques  along  with  modeled  and  approximated 
disturbance  torques.  The  relative  appendage  motion  also  causes  changes  in  the 
spacecraft’s  moments  of  inertia  used  in  the  dynamic  model.  The  kinematics  subsystem 
propagates  the  actual  spacecraft  attitude  quaternions  referenced  to  an  inertially  fixed 
coordinate  frame. 

C.  CONTROL  SYSTEM 

Primary  spacecraft  attitude  control  is  conducted  through  momentum  exchange 
with  reaction  wheels.  The  reaction  wheel  commands  are  based  on  control  laws  and  up- 
linked  feed  forward  torque  command  profile.  The  control  laws  use  errors  between  the 
measured  and  commanded  (desired)  spacecraft  angular  rates  and  attitude  quaternions. 
The  relative  momentum  generated  in  the  wheels  is  subtracted  from  the  spacecraft 
momentum  in  the  dynamics  subsystem.  Magnetic  torques  are  generated  based  on 
reaction  wheel  momentum  build  up  when  the  system  is  set  for  momentum  dumping. 
These  torques  are  directly  applied  in  the  spacecraft  dynamics  block. 

D.  ATTITUDE  DETERMINATION  SYSTEM 

The  attitude  determination  subsystem  is  of  primary  concern  in  this  thesis.  It 
incorporates  the  dynamic  modeling  concept  and  an  error  state  Kalman  filter  in  order  to 
correct  for  attitude  propagation  errors.  The  option  of  using  modeled  mechanical  gyros  to 
determine  spacecraft  angular  rate  is  also  maintained  in  order  to  allow  comparison 
simulations  to  be  conducted.  The  concept  of  using  a  dynamic  model  for  analytical  rate 
determination  will  be  referred  to  in  the  sequel  as  the  dynamic  gyro. 

The  attitude  determination  computer  simulation  uses  a  trigger  to  control  the 
bandwidth  of  discrete  calculations.  The  basic  data  flow  within  the  attitude  determination 
computer  is  shown  in  the  large  rectangular  subsystem  block  in  Figure  2.  The  equivalent 
SIMULINK  diagram  that  implements  this  subsystem  is  shown  in  Figure  4. 
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Figure  4.  SIMULINK  Attitude  Determination  Subsystem 


Random  noise  and  drift  rates  are  added  to  the  actual  angular  velocity  vector  for 
the  simulated  rate  input  of  the  mechanical  gyro  option.  Bias  corrections  from  the  Kalman 
estimator  are  applied  to  the  angular  rate  before  it  is  used  in  the  attitude  propagator.  When 
the  dynamic  gyro  option  is  simulated  the  angular  rate  is  analytically  determined  from 
measured  and  known  spacecraft  parameters  fed  into  a  discrete  dynamic  model.  Artificial 
errors  and  noise  are  applied  to  measured  and  derived  parameters  though  internal  sensor 
models.  Momentum  corrections  to  the  dynamic  gyro  are  derived  from  Kalman  filter 
updates.  The  calculated  rate  is  then  supplied  to  the  attitude  propagator. 

The  error  state  Kalman  filter  algorithm  depends  on  updates  generated  by  star 
tracker  measurements.  A  star  tracker  model  is  used  to  produce  artificial  noise  corrupted 
horizontal  and  vertical  measurements  related  to  the  star  tracker  orientation.  The  Kalman 
filter  tracks  rate  errors  for  bias  correction  of  gyro  measurements  or  momentum 
corrections  for  the  dynamic  gyro.  It  also  produces  attitude  corrections  that  are  applied  to 
the  attitude  propagator  output. 
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The  attitude  Propagator  uses  discrete  kinematics  to  convert  angular  rate  to 
quaternion  attitude.  These  estimated  parameter  are  then  used  by  the  attitude  control 
algorithm  to  complete  the  feedback  loop. 


16 


III.  DYNAMICS  AND  CONTROL 


As  three- axis  stabilized  spacecraft  become  more  technologically  advanced,  their 
operations  require  more  slewing  maneuvers  and  their  dynamic  complexities  increase. 
This  often  leads  to  multiple  rigid  or  flexible  components  that  have  independent  pointing 
and  tracking  requirements.  Complex  satellites  often  consist  of  a  primary  payload  that 
demands  strict  pointing  control  while  directional  telemetry  and  command  antennas  or 
secondary  payloads  are  controlled  independently  for  other  functions  such  as  tracking  a 
ground  station  throughout  its  maneuvers. 

Here  we  consider  the  dynamics  of  a  spacecraft  that  consists  of  a  primary  body 
with  momentum  exchange  control  devices  and  a  coupled  rigid  secondary  body  or 
appendage.  The  secondary  body  rotates  with  one  degree  of  freedom  relative  to  the 
spacecraft  about  an  axis  through  the  centers  of  mass  of  both  bodies.  Under  these 
conditions,  the  position  of  the  spacecraft  center  of  mass  remains  stationary  during 
appendage  relative  motion. 

A.  BIFOCAL  RELAY  MIRROR  SATELLITE  DYNAMICS 

The  example  spacecraft  for  which  the  equations  of  motion  are  derived  is  the 
Bifocal  Relay  Mirror  satellite  shown  in  Figure  5.  This  spacecraft  is  designed  to  receive 
laser  energy  from  a  ground  station  through  a  receive  telescope  and  to  redirect  it  through 
an  optically  coupled  transmit  telescope  to  a  different  point  on  the  ground.  The  primary 
body  is  the  transmit  telescope  and  the  rigidly  attached  spacecraft  bus  while  the  smaller 
receive  telescope  is  treated  as  the  secondary  body  or  appendage.  Weight  and  balance 
design  ensures  that  the  centers  of  mass  of  the  two  bodies  lie  close  to  the  coupled  axis  of 
rotation  so  that  the  system  center  of  mass  is  nearly  constant  during  relative  motion. 
Pointing  control  of  the  spacecraft  is  accomplished  with  reaction  wheels  while  a  drive 
motor  is  used  to  control  the  relative  angle  between  the  transmit  and  receive  telescopes. 
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There  are  three  coordinate  systems  defined  for  the  development  of  the  dynamic 
equations.  These  coordinate  frames  are  depicted  in  the  system  diagram  shown  in  figure 
x.  The  xyz  coordinate  system  is  fixed  to  the  primary  body  with  the  origin  at  its  center  of 
mass.  The  x-axis  is  oriented  parallel  to  the  axis  of  rotation  between  the  primary  and 
secondary  bodies  and  the  z  axis  is  parallel  to  the  optical  axis  of  the  telescope.  The  y  axis 
is  defined  such  that  the  xyz  coordinate  system  is  a  right-handed  mutually  orthogonal 
frame.  The  x’y’z’  coordinate  system  is  similarly  oriented  to  the  secondary  body  with  its 
origin  at  the  center  of  mass.  The  x  and  x’  axes  remain  parallel  during  appendage  motion 
while  the  angular  displacements  of  the  y’  and  z’  axes  from  the  y  and  z  axes  respectively 
are  defined  by  the  relative  rotation  angle  a.  The  equations  of  motion  are  derived  for  the 
spacecraft  body  coordinate  system,  XYZ,  which  is  parallel  to  xyz  frame.  Its  origin  is  at 
the  total  spacecraft  center  of  mass.  Unit  vectors  along  the  X,  Y,  and  Z  spacecraft  body 

axes  are  given  by  i  ,  j  and  k  respectively. 

1.  Rigid  Body  Equations  of  Motion 

The  angular  equations  of  motion  are  derived  from  the  application  of  Newton’s 
second  law  to  rotational  dynamics.  In  the  general  case  the  equation  of  motion  defined  in 
an  inertial  frame  for  a  rigid  body  about  an  arbitrary  point  P  is  given  by 

MP  =  fiP+ji("cf<ii)n  (3.1) 

m 

where  MP  is  the  net  external  torque  applied  to  the  body  about  P,  Hp  is  the  angular 
momentum  of  the  body  about  P,  rc  is  the  position  of  the  center  of  mass  relative  to  P,  nc  is 

the  vector  from  the  center  of  mass  to  the  position  of  dm  in  the  body,  and  dm  is  an 
incremental  unit  of  mass  within  the  body 

If  point  P  is  coincident  with  the  body  center  of  mass,  then  nc  =  0  and  the  equation 
simplifies  to 

M=H  (3.2) 

Equation  (3.2)  applies  in  inertial  reference  frames.  To  extend  it  to  body  coordinates 
where  it  can  be  employed,  we  need  to  understand  the  concept  of  vector  derivatives  in  a 
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rotating  coordinate  frame.  If  the  body  frame  rotates  relative  to  an  inertial  reference  frame 
with  angular  velocity  u  then  the  derivative  of  any  vector  A  in  inertial  coordinates  can 
be  related  to  the  derivative  in  the  rotating  body  coordinates  by 


A 


I 


+  u  xA 

B 


(3.3) 


Applying  this  relation  to  the  angular  momentum,  the  equation  of  motion  for  a  rigid  body 
in  rotating  body  coordinates  becomes 

M  =  H+ uxH  (3.4) 


2.  Multi- Body  Equations  of  Motion 

For  the  two  body  Bifocal  Relay  Mirror  system  with  reaction  wheels  shown  in 
figure  1,  the  total  system  angular  momentum  in  body  XYZ  coordinates  can  be  written 

H  =  HS_+Hrel+Hw  (3.5) 

where  H  is  the  total  system  angular  momentum  of  the  spacecraft,  Hs_  is  the  angular 
momentum  of  the  system  due  to  the  rotation  of  the  body  coordinate  frame,  neglecting  the 
contribution  due  to  relative  motion  of  the  receive  telescope  and  reaction  wheels,  Hre]  is 

the  angular  momentum  due  to  the  relative  motion  of  the  secondary  body,  and  Hw  is  the 
angular  momentum  due  to  the  relative  motion  of  the  reaction  wheels 

Since  the  Bifocal  Relay  Mirror  spacecraft  is  approximated  as  a  system  of  rigid 
bodies,  we  can  substitute  the  total  system  angular  momentum,  H ,  into  Equation  (3.4)  to 
get  the  multi-body  spacecraft  equation  of  motion 

M  =Ss.  +Srel  +  Sw  +  ux(Hs.  +  Bel  +  !!„  )  (3.6) 

where  M  is  the  net  external  torque  applied  to  the  spacecraft  about  its  center  of  mass 
including  all  control  and  disturbance  torques. 
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3. 


Moments  of  Inertia 


To  determine  these  angular  momentums  we  need  to  consider  the  moments  of 
inertia  of  the  spacecraft  and  it  components.  We  define  IT  to  be  the  inertia  matrix  for  the 
primary  body  (transmit  telescope  and  bus)  about  its  center  of  mass  in  the  xyz  coordinate 
frame. 
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(3.7) 


The  vector  from  the  center  of  mass  of  the  primary  body  to  the  spacecraft  center  of  mass  is 
given  by 

rT  =  x  T;  +  y,  j  +  z.yk  (3.8) 


The  moment  of  inertia  matrix  for  the  receive  telescope,  Ir,  about  its  center  of  mass  given 
in  the  x’y’z’  coordinate  frame  is  given  by 
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(3.9) 


The  vector  from  the  center  of  mass  of  the  receive  telescope  to  the  spacecraft  center  of 
mass  is  given  by 

rR  =  x  Ri  +yRj+zRk  (3.10) 


To  align  the  receive  telescope  inertia  matrix  with  the  spacecraft  body  coordinate 
frame  a  transformation  matrix  is  applied.  This  time  varying  matrix  is  the  direction  cosine 
matrix  that  defines  a  single  rotation.  The  axis  of  rotation  is  parallel  to  the  body  X  axis 
and  is  of  magnitude  a.  The  transformation  matrix  is  given  by 


1  0  0 


T  = 


0  cos(a)  sin(a) 
0  -sin(a)  cos(a) 


(3.11) 
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The  SIMULINK  subsystem  used  to  generate  this  x-axis  rotation  matrix  given  an  input 
angle  is  shown  in  Figure  6. 


Figure  6.  SIMULINK  Subsystem  Diagram:  X-Axis  Rotation  Matrix 


To  obtain  the  total  system  inertia  matrix  about  the  spacecraft  center  of  mass  in  the 
body  coordinate  frame  XYZ  the  parallel  axis  theorem  is  applied  to  the  inertia  matricies  of 
each  body.  Additionally  the  rotation  matrix  is  applied  to  the  receive  telescope  inertia 
matrix  to  align  it  with  the  body  coordinates.  If  mr  and  rp<  are  the  masses  of  the  transmit 
telescope  and  receive  telescope  respectively  the  total  system  inertia  matrix  is  given  by 
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(3.12) 


The  calculation  of  the  spacecraft  system  inertia  matrix  is  implemented  in  the  SIMULINK 
subsystem  shown  in  Figure  7. 
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Figure  7.  SIMULINK  Subsystem  Diagram:  Spacecraft  Moment  of  Inertia  Matrix 
a.  Rate  of  Change  of  Spacecraft  Inertia  Matrix 

All  components  of  the  system  moment  of  inertia  matrix  are  constant  except 
TxtIrTx  which  varies  with  the  relative  rotation  angle  a  .  Therefore  the  rate  of  change  of 
the  system  inertia  matrix  is  given  by 

i  =  txTIRTx+TxTIRTx  (3.13) 
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shown  that  I  as  a  function  of  a  and  a  is  given  by 
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(3.14) 


The  SIMULINK  subsystem  diagram  that  calculates  the  rate  of  change  of  the  spacecraft 
inertia  matrix  based  on  appendage  relative  angular  orientation  and  rate  is  shown  in  Figure 
8. 
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Figure  8.  SIMULINK  Subsystem  Diagram:  Rate  of  Change  of  Spacecraft  Inertia 

Matrix 

4.  Angular  Velocity 

The  spacecraft  angular  velocity  is  defined  by  the  angular  velocity  of  the  primary 
body  with  respect  to  the  inertial  reference  frame.  This  angular  velocity  vector  is  given  by 

Ux' 

uT=  uy  (3.15) 

_uz_ 

The  relative  angular  velocity  between  the  primary  and  secondary  bodies  due  to  the 
rotation  about  the  body  x  axis  is  given  by 
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(3.16) 


The  angular  velocity  of  the  receive  telescope  is  equal  to  the  spacecraft  angular  velocity 
plus  the  relative  velocity 


"  U  x  +Ure, 


u  x  +  a 
tly 

u. 


(3.17) 


5.  Angular  Momentum 

The  components  of  the  total  spacecraft  angular  momentums  given  in  Equation 
(3.5)  can  now  be  defined.  The  system  angular  momentum  neglecting  the  momentum  due 
to  the  relative  motion  of  the  receive  telescope  and  the  reaction  wheels  is  given  by 

Hs_=r  (3.18) 


The  angular  momentum  of  the  receive  telescope  relative  to  the  spacecraft  is  given  by 


Hrel  Ir  rel' 


(3.19) 


Substituting  Equations  (3.18)  and  (3.19)  into  Equation  (3.5)  we  get  the  total  system 
angular  momentum 

H  =  I-hu  Rlre(i+“J:  (3.20) 


Substituting  this  relation  into  Equation  (3.6),  the  spacecraft  equation  of  motion  can  be 
rewritten  as 


M=i-4i  "ft  RTreq+"jfr->h(  ~Ki  R IreV+  ~  j) 


(3.21) 


6.  Solving  for  Spacecraft  Angular  Rate 

To  solve  for  the  spacecraft  angular  rate,  u  ,  we  can  isolate  u  in  Equation  (3.21) 
and  perform  the  integration  using  a  computer  solver.  This  however  places  an 
unnecessary  burden  on  the  processor  to  continuously  calculate  the  time  derivative  of  the 

spacecraft  moment  of  inertia.  A  simpler  method  is  accomplished  by  first  solving  for  Hs_ . 
The  spacecraft  equation  of  motion  is  rewritten  as 
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M=Hs_+Ir-J.+ 


(3.22) 


which  can  be  solved  for  Hs  . 


Hs  =  M  -IR\,  -  n  +  tk “U  (3.23) 

After  the  integration, 

+  MH  H<-L)d(-  {T^,-  a  (3.24) 

t  t 

the  angular  rate,  u  ,  is  obtained  from 

u  =  !%_  (3.25) 

Figure  9  shows  the  SIMULINK  diagram  for  implements  the  spacecraft  dynamics  for  the 
Bifocal  Relay  Mirror  attitude  simulation.  The  integration  is  performed  using  the 
Doimand-Prince  ode5  solver. 


Figure  9.  SIMULINK  Subsystem  Diagram:  Spacecraft  Dynamics 

B.  COMMAND 

Maneuvering  spacecraft  often  require  feed  forward  command  to  maintain  precise 
tracking  requirements  throughout  their  maneuvers.  The  Bifocal  Relay  Mirror  satellite 
must  maintain  tightly  controlled  ground  tracking  by  both  the  transmit  and  receive 


26 


telescopes  in  order  to  perform  its  mission.  For  this  system  the  envisioned  feed  forward 
command  will  include  control  torques  necessary  to  maintain  a  dynamic  attitude  profile, 
which  includes  the  relative  motion  between  the  primary  and  secondary  bodies,  calculated 
in  the  absence  of  disturbances.  The  spacecraft  control  laws  will  be  used  to  correct  for 
errors  in  the  calculated  command  profile  and  external  disturbances  based  on  errors  in 
spacecraft  attitude  and  angular  rate  from  the  attitude  determination  system. 

The  estimated  spacecraft  rotational  profile  is  predetermined  prior  to  the  execution 
of  a  maneuver  and  tracking  operation.  The  profile  includes  the  spacecraft  body  attitude, 
angular  rate  and  angular  acceleration  as  well  as  the  relative  angle,  rate  and  acceleration 
between  the  transmit  and  receive  telescopes.  The  net  external  torque  required  to  maintain 
this  profile  in  the  absence  of  disturbances  is  fed  forward  to  the  control  devices.  This 
torque  can  be  calculated  from  the  spacecraft  equation  of  motion,  Equation(3.21), 
neglecting  the  reaction  wheel  control  devices 

Mc=i“c-ii  \H1  R~yn-“C><1  (3.26) 

In  this  equation  the  subscript  c  is  used  to  represent  feed  forward  command.  The 
equivalent  SIMULINK  subsystem  used  in  the  simulation  to  generate  the  feed  forward 
command  torque  is  shown  in  Figure  10. 


Figure  10.  SIMULINK  Subsystem  Diagram:  Feed  Forward  Torque  Command 
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C.  CONTROL 

Primary  attitude  control  of  the  Bifocal  Relay  Mirror  satellite  is  accomplished  with 
reaction  wheels.  Four  reaction  wheels  are  arranged  in  a  pyramid  constellation  to  achieve 
redundancy  in  the  event  of  a  single  wheel  failure.  Under  normal  operations  three  wheels 
are  operating  while  the  forth  is  off  line.  Control  torques  commands  are  calculated  in 
body  coordinates  and  then  distributed  to  the  three  operating  wheels. 

Torque  commands  are  generated  from  the  feed  forward  command  profile  plus  the 
attitude  control  laws  with  compensation  for  the  gyroscopic  torques  generated  from  the 
spinning  reaction  wheels.  The  simple  partial  plus  derivative  type  controller  is  chosen  for 
this  satellite.  The  control  laws  are  based  on  attitude  quaternion  and  angular  rate  error  as 
calculated  from  the  outputs  of  the  attitude  determination  system  where  the  quaternion 
error  is  calculated  by  Equation  (4.16)  [Ref.  9].  The  three  body  axis  wheel  control  laws 
are  given  by 

Hwi  =  KqlqE1Cft  h 

Hw2  =  Kq2qE2ai^2  u  (3.27) 

HW3  —  ^-q3(lE3u'i'  ^3  ^ 

More  exotic  control  schemes  exist  but  the  optimal  state  estimates  provided  by  the 
attitude  determination  system  are  simply  employed  by  these  control  laws. 

The  control  law  gains  are  chosen  to  minimize  steady  state  errors  attitude  errors 
while  ensuring  that  oscillations  induced  by  the  attitude  control  system  do  not  interfere 
with  on-orbit  structural  modes  and  payload  components.  To  determine  the  optimal 
control  gains,  attitude  control  simulations  are  conducted  with  a  representative 
maneuvering  profiles  and  external  disturbances.  The  state  inputs  used  by  the  controller 
during  gain  adjustment  simulations  are  perfect  attitude  and  rate  knowledge.  As  explained 
in  Chapter  I,  the  optimal  controller  can  be  determined  independently  from  the  attitude 
determination  system  since  its  outputs  are  based  on  optimal  state  estimation. 

The  reaction  wheel  control  subsystem  used  in  the  SIMULINK  attitude  simulation 
model  is  shown  in  Figure  11.  Saturation  and  rate  limiting  is  applied  to  simulate  nonlinear 
effects  in  real  reaction  wheel  control  systems.  The  net  control  torque  is  determined  and 
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then  distributed  to  the  individual  reaction  wheels  based  on  their  orientation  within  the 
constellation. 


D.  DISTURBANCE  TORQUES 

Representative  disturbance  torques  are  simulated  in  order  to  observe  spacecraft 
attitude  performance  in  a  realistic  space  environment. 


1.  Gravity  Gradient  Torque 

At  low  altitudes  the  torque  induced  by  gravity  gradient  on  spacecraft  without 
matched  body  moments  of  inertia  can  be  significant.  A  model  for  the  gravity  gradient 
torque  on  a  spacecraft  in  orbit  is  given  by 

MG=-VRoxIR0  (3.28) 

where  Ro  the  distance  from  the  spacecraft  center  of  mass  to  the  Earth’s  center,  R0  is  the 
unit  vector  in  that  direction  given  in  body  coordinates,  and  I  is  the  total  spacecraft 
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moment  of  inertia  matrix.  If  the  direction  cosine  matrix  from  the  orbit  to  body  frame, 
BC°  ,  is  known,  R()  is  given  by 


Rn  =  BC° 
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In  the  attitude  simulation  model,  the  orbital  reference  frame  direction  cosine  matrix 
components  are  propagated  using  the  spacecraft  angular  rate  relative  to  the  rotating  orbit 
frame  from  an  initial  orientation.  The  relative  orbital  rate  is  determined  from  the  inertial 
rate  minus  the  rate  of  rotation  of  the  orbit  frame.  The  SIMULINK  subsystem  for 
propagating  the  orbital  reference  is  shown  in  Figure  12. 


Figure  12.  SIMULINK  Subsystem  Diagram:  Orbital  Reference  Propagation 
Using  the  orbital  reference  coordinates  the  gravity  gradient  torque  can  then  be 

written 
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(3.30) 


Figure  13  shows  the  SIMULINK  subsystem  that  models  the  gravity  gradient  torque  based 
on  the  orbital  reference  frame  direction  cosine  matrix. 


Figure  13.  SIMULINK  Subsystem  Diagram:  Gravity  Gradient  Torque  Model 
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2.  Other  External  Disturbance  Torques 

Other  disturbance  torques  include  those  due  to  unbalanced  solar  pressure, 
magnetic  interactions,  and  aerodynamic  drag  effects.  These  disturbances  are  not  as  easily 
modeled  as  gravity  gradient  but  their  characteristics  and  magnitudes  are  important 
considerations  in  the  design  of  spacecraft  attitude  control  systems.  In  the  simulation 
model,  secular  and  slowly  varying  periodic  moments  are  introduced  in  each  body  axis  to 
account  for  these  unknown  disturbances.  Worst  case  magnitudes  are  chosen  based  on 
orbit  profile  and  spacecraft  characteristics  to  ensure  robust  control  design. 

E.  MAGNETIC  MOMENTS 

The  magnetic  moment  imparted  on  the  spacecraft  by  the  earth’s  magnetic  field  is 
dependent  upon  the  magnetic  field  strength,  B ,  and  the  spacecraft’s  magnetic  dipole 
vector,  m .  The  magnetic  moment  is  given  by 

Mm  =  mxB  (3.31) 


1.  Magnetic  Field  Model 

Highly  accurate  models  of  the  earth’s  magnetic  field  have  been  developed  but  a 
simple  dipole  model  is  adequate  for  the  purposes  of  this  simulation.  This  approximation 
assumes  a  simple  dipole  magnetic  field  tilted  11  degrees  from  the  earth’s  spin  axis.  The 

earth’s  magnetic  field  is  a  function  of  the  earth’s  unit  dipole  vector,  M  ,  the  distance  from 
the  center  of  the  earth  to  the  center  of  mass  of  the  spacecraft,  R,  and  the  unit  vector  in 

that  direction,  R  .  The  magnetic  field  is  given  by 


K’  r 

B  =—  3(mDR)R-M 
R  l 


(3.32) 


where  K  is  the  earth’s  magnetic  field  constant  equal  to  7.943  x  1015  Nnr/a2.  Using 
classical  orbital  elements  it  can  be  shown  that  the  components  of  the  magnetic  field 
vector  in  orbital  coordinates  are  given  by 
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cos(y)(cos(a)sin(z)  -sin(a)cos( )' cos(  r)  )-sin(Y)sin  (a)sin (u) 
-cos(a)cos (i)  -sin(a)sin(  )‘cos(  z) 

2sin(v)  (cos(a)sin(  z)  -sin(a)cos(  )‘cos(  z)  )-£cos(\)s  in(a)sin(zz) 

(3.33) 

where  v  is  the  true  anomaly  of  the  spacecraft,  £  is  the  magnetic  dipole  tilt  angle,  i  is  the 
orbit  inclination,  and  u  is  the  angle  of  the  magnetic  dipole  normal  to  the  line  of  ascending 
nodes.  The  SIMULINK  subsystem  that  propagates  the  earth  magnetic  field  vector  with 
the  orbital  reference  is  shown  in  Figure  14. 


Br 


K 

R1 


Figure  14.  SIMULINK  Subsystem  Diagram:  Earth  Magnetic  Field  Vector  in  Orbital 

Coordinates 


The  magnetic  field  vector  can  be  transformed  to  the  spacecraft  body  coordinates 
with  the  direction  cosine  matrix. 

Bb  =  bC°B0  (3.34) 

where  BC°  is  the  transformation  matrix  from  the  orbit  coordinate  frame  to  the  body 
frame.  Figure  15  shows  the  subsystem  that  realizes  this  transformation  using  direction 
cosine  matrix  components. 
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Figure  15.  SIMULINK  Subsystem  Diagram:  Magnetic  Field  Conversion  to  Body 

Coordinates 

2.  Magnetic  Control  Torque 

Using  magnetic  torquers  a  spacecraft  magnetic  dipole  can  be  generated  to  react 
with  the  earth’s  magnetic  field  to  produce  a  control  torque.  In  this  simulation  magnetic 
control  can  be  used  to  help  desaturate  the  reaction  wheels.  The  control  law  for  magnetic 
dumping  of  reaction  wheel  momentum  is  given  by 

M..,  =  -k(BBxHw)  (3.35) 

where  k  is  the  magnetic  torquer  gain.  Figure  16  shows  the  SIMULINK  subsystem  that 
simulates  the  magnetic  dumping  control  torque.  Saturation  is  added  to  simulate  nonlinear 
effects  in  the  torque  rods. 


Figure  16.  SIMULINK  Subsystem  Diagram:  Magnetic  Torquers  for  Momentum 

Dumping 


Magnetic  torquers  can  also  be  used  for  attitude  control.  The  control  laws  for  the 
torquers  to  produce  a  desired  control  moment,  MCD ,  on  the  spacecraft  is  given  by 


M 


cm 


BbxMcd 
-  2 

Bb 


(3.36) 
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IV.  ATTITUDE  REPRESENTATION  AND  KINEMATICS 


There  are  many  ways  of  describing  the  orientation  of  one  coordinate  system 
relative  to  another.  The  most  common  descriptors  used  in  spacecraft  attitude 
determination  include  Euler  angles,  direction  cosine  matricies  and  quaternions,  also 
known  as  Euler  parameters  [Ref.  10].  Euler  angles  provide  a  convenient  way  to  represent 
attitude  and  are  usually  he  easiest  to  visualize.  However,  singularities  arise  when  the 
relative  orientation  from  the  reference  coordinate  system  becomes  large.  Therefore, 
highly  maneuverable  spacecraft  require  other  means  of  attitude  representation.  Direction 
cosine  matricies  and  quaternions  overcome  this  problem.  Direction  cosine  matricies 
provide  the  most  convenient  way  of  transforming  vectors  between  coordinate  systems  but 
require  significantly  higher  attitude  processor  throughputs  than  quaternions.  As  an 
example,  propagating  an  attitude  matrix  with  angular  rate  data  requires  the  integration  of 
nine  elements  while  the  quaternion  has  only  four.  For  these  reasons,  the  quaternion 
representation  is  chosen  for  this  attitude  simulation  model. 

A.  QUATERNION  DEFINITION  AND  CHARACTERISTICS 

The  four  vector  quaternion  based  representation  q  e  R4  can  be  written  as 

qi 

-  q2 

q  = 

q3 

_q4_ 

or  equivalently, 

q=qli+q2j  +q3fc+q4 

where  the  quaternion  has  a  three  vector  imaginary  part 
qi" 

cp  =  q2  and  a  scalar  real  part  qR  =  q4 . 

q3 


(4.1) 


(4.2) 
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1.  Meaning  of  Quaterions 

Quaternions  represent  the  angular  orientation  of  a  body  relative  to  a  reference 
coordinate  frame  by  a  single  axis  rotation  of  magnitude  e  about  the  eigenvector  axis 
given  by  e  corresponding  to  the  +1  eigenvalue  of  the  direction  cosine  or  attitude  matrix. 
The  direction  cosine  matrix  can  be  written  in  terms  of  these  four  parameters 

1  0  0  0  -  3  e  2  s 

A=  0  1  0  cos(e)+  e3  0  -el  sin(e)  +  eeT(l  -cos(e))  (4.3)  (4.4) 

001  -2lee0 

Similarly  these  parameters  can  be  represented  by  the  four  quaternions 

"#1  r  i 

qt=  q2  =5  2 

_qfJ  [  3 

qR  =q4  =COS 

which  have  the  added  property  that 

||q||  =  qf  +  q2  +  q3  +  q4  =  1  (4-6) 

2.  Attitude  Matrix 

The  equivalent  attitude  or  direction  cosine  matrix  can  be  generated  from 
quaternions  using 

qf  - q2  - q3  +  q4  2(qiq2+q3q4)  2(q,q3-q2q4) 

A  =  2(q,q2-q3q4)  -qj +q22-q^  + q42  2(q2q3  +  q,q4)  (4.7) 

_  2(q!q3  +q2q4)  2(q2q3  -q^4)  -q2  -q22  +q3  +q2_ 

The  SIMULINK  subsystem  used  in  this  simulation  model  to  convert  a  quaternion  to  the 
equivalent  attitude  matrix  is  shown  in  Figure  17. 
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Figure  17.  SIMULINK  Subsystem  Diagram:  Attitude  Matrix  from  Quaternions 

3.  Quaternion  Multiplication 

The  multiplication  of  two  quaternions  to  form  a  third  defines  an  angular 
orientation  resulting  from  two  eigen-axis  rotations.  If  q  represents  the  transformation 
from  coordinate  frame  A  to  B  and  q  represents  the  transformation  from  frame  B  to  C 
then  the  transformation  from  frame  A  to  C  is  given  by 

q  =  qq  (4.8) 

This  multiplication  can  be  implemented  several  ways.  The  two  quaternions  can 
be  multiplied  directly  using  the  quaternion  format  given  in  Equation  (4.2) 

(qp  +  q2j  +  q  3k  +  q4 )  =  (q \i  +  q2  j+q3k  +  qA  )(q'i  +  qj+q'-jc  +  q')  (4.9) 


using  the  equalities 

i2  =j 2  =k2  =  — 1 
ij  =  -ji  =  k 
jk  =  ~kj  =  i 
ki  =  -ik  =  j 


(4.10) 


Quaternion  multiplication  can  also  be  performed  by  treating  the  imaginary  and  real  parts 
separately 
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q.  =  q'Rqr+  q^q.  +qi  xq" 
qR  =qRqR-qrqI 


In  matrix  form  quaternion  multiplication  is  given  by 


q  4 

-q3 

q2 

i - 

vd* 

- 1 

hO 
^  * 

q3 

q4 

-qi 

q2 

// 

q2 

-q2 

q'i 

q4 

q3 

// 

q3 

-qi 

-q2 

-q3 

q4_ 

Lq4 

(4.11) 


(4.12) 


This  multiplication  method  is  used  in  the  attitude  simulation.  The  SIMULINK  subsystem 
is  shown  in  Figure  18. 


Figure  18.  SIMULINK  Subsystem  Diagram:  Quaternion  Multiplication 

4.  Coordinate  Rotations 


Sequential  axis  rotations  such  as  those  that  define  the  three  Euler  angle 
representation  can  be  realized  by  the  successive  quaternion  products.  Additionally,  if 
there  are  three  small  simultaneous  rotations  e, ,  cn ,  and  e3  about  the  coordinate  axes  x,  y, 
and  z  respectively  the  resulting  quaternion  is  determined  as  follows: 

Let  0  =  (e2  +  e2  +e2)1/2  and  9=e- 

0 


,  then 


4i  = 


4i 

q2 

q3 


=  0  sin 


(  0i 
v  ^  ) 


and  qR  =q4  =cos 


v  2  J 


(4.13) 
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5.  Quaternion  Inverse  and  Identity 

The  inverse  of  a  quaternion  is  its  complex  conjugate  and  is  analogous  to  the 
transpose  of  a  direction  cosine  matrix.  A  quaternion  is  conjugated  by  reversing  the  sign 
on  the  vector  part. 


qi 

-qi" 

q2 

then  q  = 

-q2 

q3 

-q3 

_q4_ 

i 

1 _ 

The  identity  quaternion  is  found  by  multiplying  any  quaternion  by  its  conjugate  as 
shown  below. 


qq*=qq  = 


0 

0 

0 

1 


(4.14) 


6.  Quaternion  Error 


A  difference  quaternion  can  be  defined  between  to  orientations  that  are  referenced 
to  the  same  coordinate  frame.  If  c\  represents  the  transformation  from  coordinate  frame 
A  to  B  and  l\  represents  the  transformation  from  frame  A  to  C  then  the  transformation 
from  frame  B  to  C  is  given  by 


q=q  q 


qi 

v  rn 

a" 

-q2 

-qi 

- 1 

& 

-  ^ 
_ i 

-qi 

q4 

qi 

-qi 

// 

q2 

qi 

-q. 

q'4 

-qi 

// 

q3 

_qi 

q2 

qi 

qi. 

Lq+J 

(4.15) 


This  can  be  used  to  calculate  the  error  quaternion  where  the  target  orientation  is 
given  by  q  and  the  actual  spacecraft  orientation  is  q  . 


—  — 

qe=q  q 


(4.16) 
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The  SIMULINK  subsystem  diagram  for  determining  the  error  quaternion  is  given  in 
Figure  19. 


Figure  19.  SIMULINK  Subsystem  Diagram:  Error  Quaternion 


The  angular  difference  in  radians  between  the  original  quaternions  is  contained  in 
the  real  part  of  the  error  quaternion. 

ee=2cos-‘(q4e)  (4.17) 

7.  Vector  Transformations 

Transforming  vectors  between  coordinate  systems  requires  two  quaternion 
multiplications  where  the  vector  is  treated  as  quaternion  with  a  real  part  of  zero.  If  q 
represents  the  orientation  of  coordinate  frame  B  with  respect  to  reference  frame  A  then  a 
vector  given  in  coordinates  of  the  reference  frame,  vA ,  can  be  transformed  to  the 
coordinates  of  frame  B  by 

vB=q*vAq  (4.18) 

B.  QUATERNION  KINEMATICS 

Actual  spacecraft  motion  is  simulated  using  the  continuous  quaternion  differential 
equations.  In  the  attitude  determination  computer  simulation  the  attitude  propagator  uses 
the  discrete  kinematic  equations. 
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1. 


Continuous  Kinematics 


The  differential  equation  for  the  quaternions  of  a  rotating  coordinate  system  can 
be  found  by  differentiating  this  equation  with  respect  to  a  fixed  reference  frame  A.  If  the 
rotational  rates  for  coordinate  frame  B  are  given  by 


u  = 


,  then  d^A 


dt 


dvn 


dt 


■0  and  ^Vb 


dt 


m  xv 


B 


The  inertial  derivative  of  equation  (4.18)  becomes 


^  _  dq*  _  _  _*_  dq 

uxvb=— vAq+q  vA— 
dt  dt 


(4.19) 


Using  the  properties  of  quaternions,  substituting  the  original  equation  for  vB  and 
noting  that  the  vector  vA  is  arbitrary  it  can  be  shown  that 


dq  1  _ 
—  =  —  qu 
dt  2 


or  in  matrix  differential  form 


dq 

dt 


<ii 

qi 

q2 

q2 

_  i 

q3 

2 

q3 

2 

_q4_ 

.04. 

3u  9 

0  Uj 

-U  j  0 

-u2  -u3 


u3 

0 


0i 

q2 

q3 

q4 


(4.20) 


(4.21) 


where  A  ( u ) 


--S-(u )  u 
-uT  0 


and  -S(u)is  the  skew  symmetric  matrix  associated  with  the 


vector  u  .  The  kinematics  subsystem  diagram  used  in  the  SIMULINK  model  is  shown  in 
Figure  20.  The  integration  is  performed  using  the  Dormand-Prince  ode5  solver. 
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Figure  20.  SIMULINK  Subsystem  Diagram:  Continuous  Kinematics 

2.  Discrete  Kinematics 


The  solution  to  the  continuous  differential  equation  when  u  is  held  constant  is 


q=e 


js,u)t 


q0: 


cos^  )i4x4+s 


-u 

v2ey 


sim 


*  ) 


(  ) 


(4.22) 


where  e 


ut 


.  This  solution  leads  to  the  discrete  implementation  of  the  quaternion 


kinematic  equations. 


The  attitude  determination  system  uses  angular  rate  information  to  kinematically 
propagate  the  spacecraft  attitude  quaternion  in  discrete  time  steps.  The  angular  rate 
vector  u  e  R3  is  integrated  with  time  step  At  to  produce  the  incremental  angle  vector 
Ae  e  R3 .  For  small  time  steps  Ae  approximates  an  eigen-axis  rotation  in  the  current 
body  coordinate  frame  so  it  can  be  related  to  the  change  in  the  attitude  quaternion  by 


Aq  =  sin 


f  Ae 

3 

f  Ae,/  +  At2j  +  At3k 

+COS 

> 

l  2A5 

l 

2  J 

) 

(4.23) 


The  SIMULINK  subsystem  used  in  the  discrete  attitude  propagator  to  obtain  the 
incremental  quaternion  step  is  shown  in  Figure  21. 
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Figure  21.  SIMULINK  Subsystem  Diagram:  Discrete  Attitude  Propagator 


The  new  attitude  quaternion  is  simply  determined  by 

qnew=qold(  “A  (4.24) 

where  the  equation  implies  quaternion  multiplication. 
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V.  KALMAN  FILTER 


The  Kalman  filter  provides  a  non-deterministic  means  of  estimating  a  system  state 
vector  using  a  state  space  mathematical  plant  model  and  sensor  measurement  data  related 
to  some  subset  of  the  state  variables.  It  is  a  stochastic  optimal  estimator  designed  to 
minimize  the  weighted  mean  square  error  in  the  state  estimate.  The  Kalman  gain  matrix 
determines  the  weighting  based  on  the  relative  confidence  between  the  past  state  estimate 
propagated  to  the  current  time  and  the  current  partial  measurement  of  state  variables. 
The  error  covariance  matrix,  the  second  statistical  moment  of  the  state  vector,  tracks  the 
confidence  in  the  state  estimate  while  a  measurement  error  covariance  matrix  relates  the 
confidence  in  the  measurement. 

A.  RECURSIVE  DISCRETE  KALMAN  FILTER 

The  Kalman  filter  can  be  implemented  recursively  since  all  of  the  information 
from  past  state  measurements  is  encapsulated  in  the  previous  state  estimate  and  error 
covariance  matrix  with  are  both  tracked.  During  state  propagation  the  error  covariance  is 
updated  to  reflect  additional  error  added  by  imperfections  in  the  plant  model.  A  recursive 
discrete  Kalman  filter  is  used  in  the  proposed  attitude  determination  scheme.  Figue  22 
illustrates  the  recursive  nature  of  the  Kalman  filter. 
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Figure  22.  Discrete  Kalman  Filter  Loop 


The  standard  equations  for  a  recursive  discrete  Kalman  filter  are  summarized  in 


Table  2  below.  For  full  background  development  and  equation  derivations  see  Appendix 
A  [Ref.  1], 


System/plant  Model 

x9  =x  k-i  kY+  ti 
wk  ~  N(0,Qk) 

Measurement  Model 

zk  =Hkxk  +vk 

vk~N(0,Rk) 

Initial  Conditions 

X0=E[x0] 

P0=E[(x0-x0)(x0-k0)T] 

Assumptions 
(uncorrelated  errors) 

E[wkVj]  =0 

for  all  j,  k 

Prediction 

State  Estimate  Extrapolation 

Error  Covariance  Extrapolation 

x9  =x  k-i"k-i 

PP=Pk-Pk-l  k-Q+  k-1 

Correction 

Kalman  Gain  Matrix 

State  Estimate  Update 

Error  Covariance  Update 

Kk=PkH^(HkPkH[+Rk)1 

Xk=Xk  +  Kk(Zk-HkXk) 

Pk+  =  (I-KkHk)Pk 

Table  2.  Discrete  Kalman  Filter  Equations 
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B.  ERROR  STATE  EXTENDED  KALMAN  FILTER  IMPLEMENTATION 


A  linear  error  state  extended  discrete  Kalman  filter  is  implemented  in  the 
simulated  attitude  determination  model  to  estimate  spacecraft  attitude  and  angular  rate. 
The  nonlinear  attitude  propagation  is  performed  discretely  outside  of  the  Kalman  filter 
according  to  Equations  (4.23)  and  (4.24)  at  the  high  bandwidth  frequency  of  the  attitude 
processor.  The  angular  rate  estimate  used  in  the  kinematic  model  is  provided  by  the 
spacecraft  gyros  or  pseudo  gyro  rate  calculations.  The  Kalman  filter  is  designed  to 
provide  bias  corrections  to  the  gyro  outputs  for  attitude  propagation  [Ref.  11].  If  the 
pseudo  gyro  rate  is  used,  the  bias  error  is  treated  as  a  rate  error  from  which  a  correction  to 
the  system  angular  momentum  is  determined.  Measurement  updates  are  provided  by  star 
trackers.  The  difference  in  the  measured  and  predicted  star  vector  is  related  to  the 
attitude  error  and  used  to  provide  corrections  to  the  state  estimate. 


1.  State  Variables 

The  Kalman  filter  used  in  this  model  estimates  six  state  variables:  a  three  vector 

of  attitude  errors,  ee  R  ,  and  a  three  vector  of  gyro  bias  errors,  b  el?1.  The  total  state 
vector  is  given  by 


x= 


(5.1) 


The  attitude  error,  e  ,  represents  the  deviation  in  the  spacecraft  attitude  relative  to 
the  inertial  reference  frame  given  by  a  vector  of  three  simultaneous  rotations.  The  bias 

error,  b ,  represents  the  change  in  bias  of  the  angular  rate  data  in  the  spacecraft  body 
coordinate  frame. 


2.  Attitude  Propagation  Error  Correction  Methods 

There  are  two  ways  to  implement  the  Kalman  filter  corrections  to  the  attitude 
propagator.  In  the  first  method,  the  attitude  propagator  is  fed  with  raw  gyro  rate 
information  and  the  Kalman  filter  maintains  the  total  gyro  bias  tracking  with  the  state 

bias  error,  b  .  In  this  method,  the  attitude  quaternion  must  be  corrected  at  each  time  step 
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by  the  filter.  The  attitude  error  vector,  e ,  is  converted  to  an  incremental  quaternion 
rotation,  Aq  and  applied  to  the  propagated  attitude  quaternion  as  shown  in  Chapter  IV 
with  Equations  (4.23)  and  (4.24)  or,  using  the  small  angle  approximation  Aq  is  given  by 


Aq=—i  +  —j  +  —k  + 1  (5.2) 

2  2  2 

Then  by  quaternion  multiplication  the  propagated  attitude  quaternion  is  updated 

qnew=(Aq)q0id  (5.3) 

The  alternative  method  of  applying  the  filter  correction  to  the  propagator  is  to  use 
the  bias  error  state,  b ,  to  correct  the  unknown  gyro  bias,  u  b  each  time  a  measurement 
update  is  obtained 


U  bnew  —  U  bold  +  ^ 


(5.4) 


For  this  method,  the  total  unknown  gyro  bias  is  tracked  separately  from  the  filter  and 
added  to  the  gyro  angular  rate  before  it  is  fed  to  the  attitude  propagator. 

Ae  =  (ug+ub)At  (5.5) 

The  attitude  error  state  correction  is  then  only  necessary  at  measurement  updates.  Using 
this  method  allows  the  Kalman  filter  error  vector  x  to  be  reset  after  each  bias 
measurement  update.  Since  the  filter  approximates  nonlinear  errors  with  a  linear  model, 
keeping  the  errors  as  close  to  zero  as  possible  improves  the  estimate.  This  Kalman  filter 
correction  method  is  used  in  the  simulated  attitude  determination  model. 


3.  Plant  Model 


A  nominal  plant  model  is  chosen  for  this  Kalman  filter  which  assumes  a  constant 
rate  bias  between  filter  updates  and  uncoupled  (linear)  attitude  errors  between 
propagation  steps.  It  is  derived  from  the  continuous  equation  of  state 
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where  A  is  the  estimate  of  the  spacecraft  attitude  or  direction  cosine  matrix.  This  linear 
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state  equation  includes  coupled  rotational  effects  through  the  time  variant  estimated 
attitude  matrix  which  changes  dynamically. 


4.  State  Transition  Matrix 


For  the  discrete  Kalman  filter  implementation  of  this  model  the  state  transition 
matrix,  ok ,  is  calculated  for  each  time  step  using  the  matrix  exponential  of  f(t)  from  the 
continuous  state  equation 


°k  =  e/(t)At  = 


'3X3 


J3X3 


AN 


3X3 


(5.7) 


where  A,1  is  the  transpose  of  the  estimated  attitude  matrix  at  the  current  discrete  time,  tk 
and  At  is  the  interval  between  steps.  The  SIMULINK  subsystem  used  to  calculate  the 
state  transition  matrix  is  shown  in  Figure  23.  The  estimated  attitude  matrix  is  calculated 
from  the  propagated  attitude  quaternion  using  Equation  (4.7). 
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Figure  23.  SIMULINK  Subsystem  Diagram:  State  Transition  Matrix 


5.  Kalman  Filter  Prediction  Equations 


The  state  transition  matrix  allows  the  state  estimate,  x ,  and  associated  error 
covariance  matrix,  P,  to  be  propagated  forward  in  the  discrete  Kalman  filter  prediction 
step 


Ai=  A  (5.8) 

Pk+i=  k°kP|P  +  (5-9) 

The  subscript  k+1  indicates  next  discrete  time  step,  f+i,  and  the  superscript  -  indicates 
predicted  future  value  based  only  on  information  up  to  the  current  time  step  tk.  The 
covariance  matrix  relates  the  confidence  in  the  associated  state  estimate.  A  larger  P 
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indicates  less  confidence  which  means  that  the  attitude  determination  system  is  assumed 
to  have  larger  errors.  Measurement  updates  taken  when  P  is  large  will  greater  impact  on 
the  estimate  than  those  taken  when  the  confidence  is  high.  The  Kalman  filter  prediction 
step  is  implemented  in  the  SIMULINK  model  as  shown  in  Figure  24. 


Figure  24.  SIMULINK  Subsystem  Diagram:  Kalman  Filter  Prediction  Step 

6.  Plant  Noise  Covariance 


The  plant  noise  covariance,  Q,  is  a  positive  definite  matrix  that  characterizes  the 
plant  error  accumulated  through  the  time  step.  At ,  assuming  that  it  can  be  modeled  as 
Gaussian  white  noise  meaning  normally  distributed  with  zero  mean.  Q  is  normally  taken 
to  be  a  diagonal  matrix  meaning  that  there  is  no  known  correlation  between  the  errors  of 
the  six  state  variables.  In  terms  of  the  variances  of  the  state  variables  the  plant  noise 
covariance  matrix  is  given  by 


6^  0  0  0  0  0 

06  220  0  0  0 
0  0  4  0  0  0 
0  0  0  ^00 

0  0  0  0  4  0 

0  0  0  0  0  4 


(5.10) 


If  the  time  step  is  varied  or  the  plant  error  is  known  to  be  significantly  changed 
due  to  operating  mode,  then  Q  should  be  modified  to  reflect  that  change.  Decreasing  Q 
effectively  asserts  that  the  plant  generates  less  error  and  the  predicted  state  estimate  will 
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be  weighted  higher  relative  to  the  measurement  updates.  On  the  other  hand,  if  it  is 
increased  less  confidence  is  placed  on  the  estimate  predicted  by  the  model  and  more 
significant  corrections  will  be  applied  during  measurement  updates.  In  this  simulation,  Q 
is  taken  to  be  constant. 

7.  Kalman  Filter  Initialization 

In  order  to  initiate  the  Kalman  filter  algorithm,  an  initial  state  estimate,  x;  and  its 
associated  covariance,  P„  must  be  chosen.  Least  squares  batch  processing  can  be 
performed  on  two  or  more  star  tracker  measurements  prior  to  initiating  the  filter  or  an 
educated  guess  can  be  used.  Since  the  state  represents  errors  from  truth,  the  selected 
initial  state  estimate  is 
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constant,  n,  is  chosen  to  reflect  confidence  in  the  estimate  at  initialization.  The  attitude 
errors  generated  are  normally  about  an  order  of  magnitude  up  from  the  rate  bias  errors  so 
the  variances  are  weighted  higher. 

8.  Sensor  Measurement  Update 

The  state  estimate  and  the  covariance  matrix  propagate  forward  without 
correction  through  each  time  step  until  a  measurement  update  is  produced.  The  state 
estimate  is  used  to  kinematically  propagate  the  attitude  quaternion  while  the  covariance 
matrix  builds  up  due  to  the  added  error  through  each  step. 

a.  Measurement  Vector 

If  an  attitude  measurement  update  is  produced  at  the  next  time  step,  k+1, 
by  one  or  more  star  trackers  the  Kalman  filter  correction  step  is  applied.  The  star  trackers 
produce  horizontal  and  vertical  outputs  (H,V)  corresponding  to  the  position  of  the  star  on 
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the  detector  array.  These  outputs  are  internally  compensated  for  with  software  and 
calibrated  to  form  a  measurement  vector  in  tracker  coordinates 


s 
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H 

V 

1 


with  the  normalization 


(5.11) 


which  is  the  star  tracker  measurement  vector.  For  simulation  modeling  purposes,  noise  is 
added  to  the  H  and  V  components  to  create  an  artificial  star  tracker  measurement  vector. 
The  measurement  vector  is  generated  in  the  SIMULINK  model  as  shown  in  Figure  25. 


Figure  25.  SIMULINK  Subsystem  Diagram:  Star  Tracker  Measurement 


At  every  integration  cycle,  At ,  the  attitude  estimate  is  available  from  the 
Kalman  filter  prediction  step.  If  a  star  tracker  observation  has  occurred  during  that  time 
step,  the  measurement  vector  must  be  propagated  ahead  to  correspond  to  the  current 
discrete  time  of  the  attitude  computer.  Otherwise,  the  state  vector  and  the  state  transition 
and  covariance  matricies  must  be  interpolated  to  accomplish  the  filter  update.  Tracker 
processing  latencies  and  transport  delays  must  also  be  compensated  for  in  the  time 
difference.  For  simplicity  in  this  model,  it  is  assumed  that  the  all  observations  occur  at  an 
exact  discrete  time  step  of  the  integration  cycle.  To  produce  star  measurements  for  the 
model,  an  artificial  star  tracker  reference  is  chosen  with  representative  noise  applied  to 
the  H  and  V  outputs. 

b.  Predicted  Measurement 

The  predicted  vector  in  star  tracker  coordinates  is  needed  to  determine  the 
measurement  residual  for  the  update.  This  is  generated  using  the  known  position  of  the 
star  in  inertial  space.  The  detected  star  goes  through  the  identification  process  and  gets 
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compensated  for  aberration  to  yield  a  unit  vector,  3,  in  inertial  coordinates.  Using  the 
estimate  for  the  inertial  to  body  attitude  matrix  generated  from  the  filter  state  prediction, 

Ak+1  ,  and  the  calibrated  body  to  star  tracker  transformation  matrix,  T,  a  predicted  vector 
is  generated  in  tracker  coordinates. 

sp=TAk+1s,  (5.12) 


The  inertial  star  vectors  for  the  model  are  generated  by  applying  the 
transpose  of  the  body  to  star  tracker  and  true  attitude  matricies  to  the  same  star  tracker 
reference  used  to  generate  the  measurement. 


Si=Ak+1T 
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(5.13) 


Figure  26  illustrates  the  SIMULINK  subsystem  used  to  produce  a  simulated  inertial  star 
reference  vector. 


Figure  26.  SIMULINK  Subsystem  Diagram:  Inertial  Star  Vector 
c.  Measurement  Residual 

The  measurement  residual,  zk+i,  is  formed  by  subtracting  the  predicted 
from  measured  star  vector  and  considering  only  the  first  two  components. 

zk+i=E(sm-sp)  (5.14) 


where  E= 


1  0  0 
0  1  0 


The  measurement  residual  is  determined  in  the  SIMULINK  model 


as  shown  in  Figure  27. 
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Figure  27.  SIMULINK  Subsystem  Diagram:  Measurement  Residual 

9.  Observation  Matrix 


Since  the  measurement  residual  is  simply  the  difference  between  unit  vectors,  it  is 
easily  related  to  the  state  attitude  errors.  The  observation  or  feedback  sensitivity  matrix 
defined  by  the  relationship  to  the  state  variables 
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k+r 


ETAt„S(s,)  0 


2x3 


(5.15) 


where  -Sfs, )  indicates  the  skew  symmetric  matrix  associated  with  the  inertial  star 
reference  vector.  Figure  28  shows  the  subsystem  diagram  used  in  the  SIMULINK  model 
to  produce  the  observation  matrix. 


Figure  28.  SIMULINK  Subsystem  Diagram:  Observation  Matrix 

10.  Kalman  Gain 


The  Kalman  gain  can  then  be  calculated  for  the  correction  step  at  k+1  from  the 
standard  discrete  filter  equation 

Kk+1=Pk+lHL(Rk+1  +Hk+1Pk+1  Hk+I  )4  (5.16) 

where  Re  Rlxl  is  the  measurement  noise  covariance  matrix  associated  with  the  assumed 
Gaussian  white  noise  in  the  H  and  V  outputs  of  the  star  trackers.  Since  these  errors  are 
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usually  similar  and  uncorrelated,  the  R  matrix  is  normally  taken  to  be  a  multiple  of  the 
identity. 


R= 


0 

2 

V 


Noisier  or  less  accurate  star  trackers  will  generally  have  a  larger  R  matrix  which 
decreases  the  Kalman  gain  and  thus  provides  less  of  a  correction  to  the  state  and 
covariance  matricies.  In  reality,  the  measurement  error  may  not  be  constant  for  each 
correction  step  even  when  the  same  star  tracker  is  used.  Higher  intensity  stars,  although 
easier  to  detect  and  identify,  produce  slightly  larger  noise.  Also  stars  detected  toward  the 
field  of  view  limits  of  the  tracker  usually  have  larger  errors  due  to  distortions  than  those 
detected  near  the  star  tracker  optical  axis. 

The  Kalman  gain  SIMULINK  subsystem  is  shown  in  Figure  29. 


Figure  29.  SIMULINK  Subsystem  Diagram:  Kalman  Gain 

11.  Kalman  Filter  Correction  Equations 

With  the  Kalman  gain,  the  state  and  covariance  matricies  can  be  corrected  with 
the  measurement  update 

A^=-Kk+1zk+1  (5.17) 

xk+J=%+  *  (5-18) 

and 


Pk+i=(/6X6-Kk+1Hk+1)Pk+1  (5.19) 

The  covariance  matrix  update  is  accomplished  in  the  simulation  model  as  shown  in 
Figure  30. 
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Figure  30.  SIMULINK  Subsystem  Diagram:  Covariance  Updeate 

In  general,  the  estimated  state  drifts  away  from  the  true  state  as  it  is  propagated 
forward  with  the  non-ideal  plant  model  and  is  corrected  back  toward  truth  as 
measurement  updates  occur.  Corrections  to  the  attitude  components  associated  with  the 
body  axes  closely  aligned  with  the  star  tracker  optical  axis  will  be  corrected  much  less 
than  those  that  are  aligned  perpendicular.  The  spacecraft  attitude  error  continues  to  grow 
during  periods  where  there  are  no  cataloged  stars  in  the  sensor  field  of  view  or  when  the 
tracker  information  is  unavailable.  Additionally,  when  only  one  tracker  provides  updates 
for  extended  periods  the  angular  orientation  about  that  star  tracker’s  axis  in  body 
coordinates  goes  unchecked. 

C.  CONTROLLER  DESIGN  IMPLICATIONS 

Since  the  Kalman  filter  is  an  optimal  least  squares  estimator,  the  development  of 
an  optimal  controller  can  be  accomplished  independently.  Therefore,  the  attitude 
determination  output  should  not  affect  the  controller  design.  In  this  model,  quaternion 
and  rate  error  control  laws  are  used  as  well  as  feed  forward  torque  to  generate  reaction 
wheel  commands.  The  optimal  gains  for  the  controller  as  determined  by  simulation  with 
ideal  deterministic  attitude  knowledge  remain  optimal  with  the  attitude  determination 
based  on  Kalman  filter  state  estimates. 


56 


VI.  DYNAMIC  RATE  CALCULATION 


The  continuous  dynamic  equations  of  motion  for  the  Bifocal  Relay  Mirror 
spacecraft  are  derived  in  Chapter  III.  These  equations  produce  the  spacecraft  angular  rate 
from  external  control  and  disturbance  moments  applied  to  the  body.  A  similar  discrete 
model  can  be  applied  in  the  spacecraft  attitude  processor  software  to  produce  a  real  time 
calculated  estimate  of  the  angular  rate,  referred  to  as  the  dynamic  gyro.  This  process  is 
borrowed  from  The  Aerospace  Corporation’s  “Pseudo  Gyro”  concept  [Ref.  8].  At  high 
bandwidth  processor  execution  the  discretization  of  the  dynamics  introduces  little  error. 
The  angular  rate  generated  by  this  method  can  be  used  as  a  substitution  for  conventional 
gyroscope  outputs.  Attitude  determination  based  on  the  dynamic  gyro  can  be 
implemented  as  a  back  up  failure  mode  or  a  primary  operating  mode  to  increase  the 
expected  lifetime  of  the  satellite  gyroscopes. 

A.  DISCRETE  EQUATIONS  OF  MOTION 

The  discretized  equations  of  motion  are  derived  from 

AH  =  ^MAt  (6.1) 

where  is  the  sum  of  external  moments  applied  to  the  spacecraft  including  controls, 

modeled  disturbances  and  gyroscopic  stiffness.  This  allows  the  total  system  angular 
momentum  to  be  tracked  with 

Hk+1=Hk+AH  (6.2) 

Then  subtracting  the  relative  momentum  of  the  reaction  wheels  and  secondary  body 
produces 

HS-=H-Hw-Hffil  (6.3) 

The  calculated  spacecraft  angular  rate  is  then  given  by 

u  =  1%_  (6.4) 
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B.  MOMENTUM  CORRECTION  FROM  KALMAN  FILTER  UPDATES 

The  accuracy  of  the  dynamic  angular  rate  calculation  ultimately  depends  on  the 
tracking  of  the  system  angular  momentum.  If  uncorrected,  the  numerous  error  sources  in 
the  model  will  cause  the  angular  rate  error  to  grow  over  time.  The  error  state  Kalman 
filter  designed  for  gyro-based  attitude  determination  systems  can  be  used  to  provide  the 

necessary  model  corrections.  The  gyro  bias  error  states,  b ,  are  interpreted  as  spacecraft 
body  rate  errors.  Using  the  calculated  spacecraft  inertia  matrix,  I,  a  correction  to  the 
system  moment  of  inertia  can  be  generated  by 

AHcorr  =  Ib  (6.5) 

The  Kalman  filter  momentum  correction  is  applied  as  if  the  error  in  the  dynamic  gyro  is 
attributable  to  the  total  spacecraft  body.  The  relative  momentum  terms  from  the 
secondary  body  and  the  reaction  wheels  are  treated  as  if  they  are  without  error. 

The  SIMULINK  dynamic  gyro  subsystem  used  in  the  attitude  simulation  model  is 
illustrated  in  Figure  31. 


Figure  31.  SIMULINK  Subsystem  Diagram:  Dynamic  Gyro 


C.  INPUTS  AND  ERROR  SOURCES 

After  initial  calibration,  kinematic  plant  error  in  gyro-based  attitude  determination 
systems  is  almost  entirely  attributable  to  a  single  set  of  imperfect  gyroscope  rate  sensors. 
As  long  as  gyro  data  does  not  become  erratic  a  Kalman  estimator  based  on  a  slowly 
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changing  rate  bias  plant  model  produces  an  effective  attitude  determination  system  even 
with  relatively  noisy  rate  inputs. 

The  error  in  rate  calculations  from  dynamic  modeling,  on  the  other  hand,  is  due  to 
numerous  factors  and  is  much  harder  to  characterize.  There  are  multiple  internal  sensors 
involved  as  well  as  dynamic  modeling  simplifications.  Since  the  dynamic  calculation  is 
produced  from  total  system  momentum  tracking  any  error  in  knowledge  of  external 
torques  directly  correlates  to  rate  error.  Errors  in  system  or  component  moments  of 
inertia  have  the  same  effect.  Like  gyro  outputs,  internal  position  and  rate  sensor  data  are 
corrupted  by  measurement  and  alignment  errors.  These  data  from  all  moving  appendages 
and  momentum  exchange  devices  are  critical  to  the  accuracy  of  the  rate  calculation. 
Satellites  not  designed  to  use  dynamic  rate  calculations  for  attitude  determination  are 
usually  not  equipped  with  appendage  relative  rate  measurement  sensors.  These  data  must 
either  be  derived  adding  more  error  to  the  calculation  or  substituted  with  commanded 
rates.  It  is  important  that  all  known  biases  be  removed  from  sensor  data  and  calculated 
input  errors  since  the  Kalman  filter  estimator  is  based  on  the  assumption  of  uncorrelated 
zero-mean  Gaussian  noise.  Even  if  all  input  parameters  were  known  exactly  the  discrete 
modeling  of  the  spacecraft  dynamics  introduces  some  error. 

1.  External  Control  and  Disturbance  Torques 

In  the  dynamic  gyro,  known  externally  applied  moments  are  integrated  in  the 
system  angular  momentum  calculations.  These  include  control  moments  other  than  those 
imparted  by  momentum  exchange  devices  as  well  as  modeled  disturbance  torques.  Since 
control  torques  are  normally  of  significant  magnitudes,  it  is  essential  that  they  be 
modeled  correctly.  Moments  from  magnetic  torque  rods  depend  on  the  imperfectly 
controlled  magnetic  dipole  and  the  earth’s  magnetic  field  strength.  The  magnetic  dipole 
is  provided  by  torquer  current  measurement.  The  earth’s  magnetic  field  must  either  be 
modeled  or  measured  with  magnetometer  but  is  not  precisely  known.  Reaction  jet 
moments  are  almost  impossible  to  model  accurately.  This  often  result  in  degraded 
pointing  and  tracking  during  firings  operations. 

The  effects  of  external  disturbance  moments  depends  on  spacecraft  configuration 

and  orbital  profile.  Those  disturbances  that  have  significant  effects  on  vehicle  dynamics 
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should  be  modeled  whenever  possible  to  increase  the  accuracy  of  the  attitude  control 
system.  These  types  of  errors  are  non- zero  mean  in  the  short  term  and  therefore  can  only 
be  corrected  for  with  sensor  measurement  updates.  An  extended  option  for  The 
Aerospace  Corporation’s  Pseudo  Gyro  includes  a  torque  bias  estimator  to  reduce  the 
effects  of  unknown  external  disturbances. 

For  the  Bifocal  Relay  Mirror  satellite,  the  gravity  gradient  torque  is  the  most 
significant  disturbance  and  can  be  modeled  as  an  input  to  increase  the  accuracy  of  the 
dynamic  rate  calculation.  The  model  is  well  understood  and  is  generated  from  the  inertia 
matrix  which  is  already  required  by  the  dynamic  gyro  and  vehicle  orientation  with 
respect  to  the  gravity  vector  which  can  be  determined  from  the  estimated  attitude  and 
ephemeris  data.  The  gravity  gradient  model  used  by  the  dynamic  gyro  is  equivalent  to 
the  subsystem  shown  in  Figure  13. 

2.  Reaction  Wheel  Relative  Momentum 

Instead  of  integrating  the  torques  produced  by  momentum  exchange  devices  their 
relative  momentum  effects  are  used  directly  in  the  dynamic  rate  calculation.  The  relative 
momentum  of  each  reaction  wheel  is  given  by  its  orientation  within  the  spacecraft,  the 
component  inertia  of  its  spinning  disk  and  the  wheel  spin  rate.  The  imperfect  sensor 
measurements  from  the  reaction  wheel  tachometers  introduce  errors  in  system 
momentum  calculation.  Relative  orientation  angles  of  reaction  wheels  are  fixed  and 
errors  can  be  corrected  through  calibration.  Orientations  of  control  moment  gyros, 
however,  are  variable.  Since  these  devices  usually  carry  much  more  momentum  small 
gimbal  resolver  errors  can  have  a  significant  impact  on  total  system  momentum 
calculations.  In  this  simulation,  time  varying  artificial  alignment  errors  are  applied  to  the 
reaction  wheel  momentum  measurements  to  observe  these  effects.  The  error  corrupted 
wheel  momentum  measurement  subsystem  implemented  in  the  SIMULINK  model  is 
shown  in  Figure  32. 
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Figure  32.  SIMULINK  Subsystem  Diagram:  Error  Corrupted  Reaction  Wheel 

Momentum  Measurement 

3.  Moment  of  Inertia  Calculations 

The  calculation  of  the  total  spacecraft  inertia  matrix  is  accomplished  by  Equation 
(3.12)  using  the  SIMULINK  subsystem  illustrated  in  Figure  7.  Fixed  component 
moments  of  inertia  and  masses  of  the  primary  and  secondary  bodies  are  assumed  to  be 
known  as  well  as  the  relative  positions  of  their  centers  of  mass  from  the  system  mass 
center.  Imperfect  knowledge  of  these  parameters  introduces  errors  in  the  angular  rate 
calculation.  Errors  generated  from  rotating  spacecraft  components  are  not  constant  in  the 
spacecraft  body  frame. 

The  inertia  matrix  also  depends  on  internal  sensor  input  from  position  encoders  or 
potentiometers  for  relative  angular  orientation  of  appendages.  The  model  of  the 
potentiometer  that  measures  the  relative  angle  of  the  receive  telescope  includes 
quantization  effects  and  additive  noise.  If  appendage  relative  motion  is  slow  or 
component  moments  of  inertia  are  small,  it  may  not  be  necessary  to  update  the  system 
inertia  dyadic  at  the  bandwidth  of  the  attitude  processor.  A  trigger  is  added  to  the 
SIMULINK  subsystem  that  calculates  the  inertia  matrix  (Figure  7),  so  that  the  affects  of 
update  rate  can  be  evaluated. 

4.  Appendage  Relative  Momentum 

Knowledge  of  appendage  relative  momentums  has  direct  bearing  on  the  system 
momentum  and  therefore  the  dynamic  rate  calculation.  The  relative  momentum  of  the 
Bifocal  Relay  Mirror  satellite’s  secondary  body  is  calculated  from  Equation  (3.19).  It 
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depends  on  the  knowledge  of  the  fixed  secondary  body  inertia  matrix  and  the  relative 
angular  rate.  The  attitude  simulation  model  assumes  that  there  is  no  directed 
measurement  of  the  relative  angular  rate.  The  rate  must  therefore  be  derived  from 
potentiometer  measurements  of  the  relative  orientation  about  the  axis  of  rotation.  This 
sensor  may  have  a  minimum  discemable  incremental  angle  and  noise  corruption.  Also, 
since  the  rate  is  derived  from  position  measurements  it  exhibits  increased  noise  and  time 
lag.  The  simulation  diagram  that  produces  the  error  corrupted  relative  angle  and  rate  for 
the  spacecraft  secondary  body  is  shown  in  Figure  33. 


Figure  33.  SIMULINK  Subsystem  Diagram:  Appendage  Relative  Angle  and  Rate 

Measurement 

Appendages  of  significant  inertia  or  relative  rates  like  the  receive  telescope  of  the 
Bifocal  Relay  Mirror  satellite  should  be  controlled  with  smooth  gimbal  drive  motors  in 
order  to  minimize  nonlinear  relative  pointing  errors.  Ideally,  all  moveable  appendages 
would  have  an  associated  rate  sensor  for  each  axis  of  rotation. 

The  other  option  for  approximating  appendage  relative  angular  rates  is  to  use  the 
commanded  drive  input.  In  some  systems,  appendage  controllers  can  provide  a  smooth 
relative  rate  through  the  drive  motor,  which  significantly  enhances  angular  momentum 
tracking.  Drive  actuators  with  considerably  erratic  friction  effects  will  cause  errors  in  the 
attitude  control  system  during  slew  operations.  Short  term  transients  may  introduce 
significant  settling  times  for  error  corrections. 

The  subsystem  of  the  attitude  determination  simulation  that  controls  all  of  the 
internal  sensor  measurements  and  input  parameter  calculations  for  the  dynamic  gyro  is 
shown  in  Figure  34. 
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Figure  34.  SIMULINK  Subsystem  Diagram:  Dynamic  Gyro  Inputs 

D.  CALCULATED  ANGULAR  RATE  ERROR  CORRECTION 

The  error  state  extended  Kalman  filter  is  designed  to  track  rate  bias  errors  that  are 
relatively  constant  in  the  spacecraft  body  frame.  This  is  a  good  approximation  for 
properly  functioning  gyro  based  systems.  For  the  filter  to  remain  effective  when  used 
with  dynamic  gyro,  the  bias  in  the  output  rate  as  seen  in  body  coordinates  due  to  the 
various  error  sources  must  be  small  and  exhibit  a  bandwidth  below  the  measurement 
update  rate.  Since  it  is  the  case  that  the  common  error  sources  are  not  zero-mean,  their 
effects  on  the  rate  output  must  be  relatively  constant  so  that  the  momentum  correction 
supplied  by  the  Kalman  filter  applies  over  the  update  interval.  Transient  error  spiking 
can  require  multiple  star  tracker  measurement  updates  to  correct.  If  the  internal  sensors 
provide  inconsistent  measurements,  dynamic  gyro  angular  rate  will  be  degraded. 
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VII.  RESULTS 


Simulation  results  demonstrate  that  multi-body  spacecraft  attitude  control  without 
the  use  of  rate  gyroscopes  can  be  performed  using  the  attitude  and  angular  rate  estimation 
scheme  proposed  in  this  thesis.  The  performance  of  dynamic  gyro  based  attitude 
determination  system  is  compared  with  a  similar  gyro  based  system  using  the  Bifocal 
Relay  Mirror  attitude  simulation.  Simulation  input  parameters  are  varied  to  analyze  the 
effects  of  major  error  sources  on  the  dynamic  gyro  model.  Also  evaluated  are  the  effects 
of  star  tracker  accuracy  and  measurement  update  rates  on  the  attitude  determination 
system. 

A.  BASELINE  SIMULATION 

A  full  set  of  MATLAB  plots  are  presented  in  this  section  for  baseline  analysis  of 
the  dynamic  gyro  based  attitude  determination  and  control  scheme  developed  in  this 
thesis.  This  set  of  results  validates  the  potential  effectiveness  of  the  proposed  attitude 
and  angular  rate  estimating  scheme  used  for  multi-body  spacecraft  control.  It  also 
provides  a  common  reference  for  analysis  and  comparison  with  subsequent  simulation 
results.  For  other  simulations,  only  selected  plots  that  are  required  for  analysis  of  results 
will  be  shown. 

1.  Simulation  Input  Parameters 

Table  3  shows  the  inputs  parameters  that  are  held  constant  for  each  simulation  run 
used  to  obtain  results.  These  input  parameters  are  set  in  the  MATLAB  script  file  that 
calls  the  SIMULINK  attitude  simulation.  The  MATLAB  code  file  is  included  as 
Appendix  B. 


SIMULATION  PARAMETERS  (SIMULINK) 

Simulation  Time  Period 

500  sec 

Attitude  Determination  Bandwidth 

20  Hz 

SIMULINK  Solver  Method 

ode5  (Dormand-Prince) 

Solver  Fixed  Step  Size 

0.05  sec 

ORBIT  PARAMETERS  (Circular  Orbit) 

Altitude 

715  km 

65 


COMMANDED  MANEUVER  PROFILE 

Inertial  Attitude  Quaternions 

See  Figure  35 

Body  Axes  Angular  Rates 

See  Figure  36 

Secondary  Body  Relative  Angle 

See  Figure  37 

DISTURBANCE  MOMENTS 

Gravity  Gradient 

Modeled 

Secular  (Magnitude) 

le-4  Nm 

Periodic  (Magnitude,  Period) 

4e-4  Nm,  [400,500,600]  sec 

Disturbance  Effect  on  System  Angular 
Momentum 

See  Figure  38. 

MOMENTS  OF  INERTIA  MATRICES 

Primary  Body 

[2997.28025,- 

3.9331,118.2824 

-3.9331,3164.18285,1.1230 

118.2824,1.1230,881.82105] 

kgm2 

Secondary  Body 

[1721.07340,-0.0116,- 

7.8530 

-0.0116,1559.85414,- 

12.5463 

-7.8530,- 

12.5463,182.89962]  kgm2 

CENTER  OF  MASS  OFFSET  FROM  SPACECRAFT  C.M. 

Primary  Body 

[0.558354158, 

3.91788e-4, 

0.15226902]  m 

Secondary  Body 

[-1.302113918, 

-9.13673e-4, 

0.355100092]  m 

WHEEL  CONTROL  LAWS 

Quaternion  Error  Gains  (Kq) 

[3000,7000,4500] 

Angular  Rate  Error  Gains  (Kw) 

[1000,2000,1000] 

Control  Law  Delay  for  Initial 

Determination  Errors 

30  sec 

GYROSCOPE  CHARACTERISTICS 

Static  Rate  Biases 

le-4*[-l, 1.5,1]  rad/sec 

Rate  Noise  Variance 

le-8 

Acceleration  Noise  Variance  (Rate 
Random  Walk) 

le-12 

3  Gyro  Alignment 

Aligned  to  Body  Axes 

ATTITUDE  DETERMINATION  U 

VITIAI  IZATION  ERRORS 

Quaternion  Errors  (ql,q2,q3) 

[0.008,0.012,-0.008] 

Angular  Rate  Errors 

[-0.001,0.001,0.002]  rad/sec 

KALMAN  FILTER  INITIALIZATION 

State  Estimate 

[0,0, 0,0,0, 0] 
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State  Error  Covariance  Matrix 

5e2*[100, 0,0, 0,0,0 
0,100,0,0,0,0 

0,0,100,0,0,0 

0,0,0, 1,0, 0,0 

0,0, 0,0, 1,0,0 

0,0, 0,0, 0,1,0 

0,0, 0,0, 0,0,1] 

STAR  TRACKER  ALIGNMENT  TO  BODY  AXES 

Tracker  1  (x-rotation,  y-rotation) 

135  deg,  30  deg 

Tracker  2  (x-rotation,  y-rotation) 

135  deg,  -30  deg 

Tracker  3  (x-rotation,  y-rotation) 

180  deg,  0  deg 

REACTION  WHEEL  PARAMETERS  (4  WHEEL  PYRAMID 
CONSTELLATION) 

Number  of  Wheels  Operating 

3 

Constellation  Angle  to  xy-Plane 

45  deg 

Constellation  Torque  Saturation  Limit 

1  Nm 

Constellation  Torque  Rate  limits 

10  Nm/sec 

Table  3.  Simulation  Input  Parameters 


2.  Command  Attitude  Profile 

The  500  second  maneuvering  profile  chosen  for  the  Bifocal  Relay  Mirror  attitude 
simulation  analysis  is  illustrated  in  Figures  35,  36  and  37.  This  profile  resembles  the 
maneuver  required  to  maintain  transmit  and  receive  telescope  pointing  control  during  an 
overhead  operational  pass  to  conduct  laser  relay  operations.  The  majority  of  the 
maneuver  is  performed  in  the  spacecraft  pitch  axis,  q2,  as  both  telescopes  orient  to  point 
at  fixed  ground  sites.  Less  significant  motion  is  required  in  the  spacecraft  roll  and  yaw 
axes  in  order  to  ensure  that  the  relative  axis  of  rotation  of  the  receive  telescope  is 
correctly  oriented  during  the  tracking  maneuver.  Based  on  ground  site  separation 
distance  and  orbital  altitude,  the  largest  relative  angle  require  between  the  telescopes  is 
about  30  degrees  during  a  near  overhead  pass  between  the  uplink  and  downlink  ground 
sites. 
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Commanded  Quaternions 


Figure  35.  Baseline  Commanded  Attitude  Profile 


0  100  200  300  400  500 


Time  (sec) 

Figure  36.  Baseline  Commanded  Angular  Rate  Profile 
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Figure  37.  Baseline  Commanded  Relative  Angle  to  Receive  Telescope  Profile 

The  magnitude  of  the  total  spacecraft  angular  momentum  during  the  maneuvering 
profile  is  plotted  in  Figure  38.  The  spacecraft  attitude  is  completely  controlled  by 
momentum  exchange  with  the  reaction  wheels.  Therefore,  the  changes  to  the  angular 
momentum  profile  are  attributable  the  external  disturbance  moments.  The  gravity 
gradient  disturbance  is  modeled  and  the  other  disturbances  are  fixed  so  the  angular 
momentum  profile  remains  essentially  common  in  all  simulation  runs. 
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3.  Attitude  Determination  System  Performance  Results 

In  this  baseline  simulation,  updates  from  randomly  selected  star  trackers  are 
provided  to  the  dynamic  gyro  and  the  attitude  propagator  via  the  Kalman  filter  at  two 
second  intervals.  The  star  tracker  H  and  V  measurements  are  corrupted  with  noise 
variance  of  lxlO"4.  In  order  to  observe  the  performance  of  the  attitude  determination 
system  without  updates,  a  200  second  star  gap  is  simulated  starting  100  seconds  into  the 
run.  As  a  worst-case  analysis,  this  star  gap  occurs  during  the  peak  maneuvering  time  of 
the  satellite  including  the  rotation  of  the  secondary  body.  Attempt  is  made  to  tailor  the 
plant  and  measurement  error  covariance  matricies  used  in  the  Kalman  filter.  The  attitude 
determination  system  is  initiated  with  the  errors  given  in  Table  3. 

The  accuracy  of  the  angular  rate  calculation  is  entirely  dependent  upon  the  ability 
of  the  dynamic  gyro  to  track  the  total  spacecraft  angular  momentum.  The  error  in  the 
magnitude  of  the  total  system  momentum  compared  to  the  simulated  actual  momentum  is 
shown  in  Figure  39.  The  steady  state  momentum  error  is  held  within  0.07  Nms  with 
consistent  star  tracker  data  but  builds  to  0.35  Nms  after  200  seconds  without  stars.  No 
external  disturbance  torques  are  modeled  as  dynamic  gyro  inputs  in  this  simulation  run. 


Figure  39.  Baseline  Dynamic  Gyro  Angular  Momentum  Error 

After  star  tracker  measurements  are  processed,  the  Kalman  filter  provides  a 
momentum  correction  to  the  dynamic  gyro  determined  from  the  spacecraft  inertia  and  the 
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rate  bias  error  state  estimate.  The  attitude  quaternion  is  also  updated  and  the  states  are 
reset  to  zero.  Figure  40  shows  the  magnitude  of  the  error  states  determined  by  the 
Kalman  filter.  No  updates  are  provided  during  the  200  second  star  gap. 
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Figure  40.  Baseline  Kalman  Filter  Attitude  and  Rate  Bias  Errors 


Quaternion  errors  in  the  attitude  determination  system  are  plotted  in  Figure  41. 
Steady  state  attitude  errors  are  maintained  within  3xl0"4  during  periods  of  continuous  star 
coverage  but  are  increased  an  order  of  magnitude  by  the  end  of  the  200  second  star  gap. 
After  the  star  gap,  the  attitude  error  build  up  is  quickly  removed  through  measurement 
corrections. 
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Quaternion  Errors 


Figure  41 .  Baseline  Estimated  Attitude  Quaternion  Error 

Figure  42  illustrates  the  nature  of  the  angular  rate  estimation  error  produced  by 
the  dynamic  gyro  based  determination  system.  The  error  generated  in  the  roll  axis,  Wi,  is 
the  most  significant  because  it  is  aligned  with  the  axis  of  rotation  of  the  secondary  body. 
The  relative  rate  is  derived  from  imperfect  potentiometer  measurements.  The  rate  errors 
increase  during  receive  telescope  motion  due  to  the  potentiometer  quatization  effect 
which  produces  the  broken  pattern  of  noise.  The  yaw  axis  is  most  coupled  dynamically 
to  roll  axis  and  exhibits  similar  errors  at  less  magnitude.  The  build  up  of  the  bias  in  the 
rate  error  is  hard  to  perceive  among  the  noise  but  the  effects  are  evident  in  the  quaternion 
error  plot. 
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Figure  42.  Baseline  Estimated  Angular  Rate  Error 


4.  Attitude  Control  System  Performance  Results 

Attitude  control  is  accomplished  by  three  reaction  wheels  selected  from  a  pyramid 
constellation  of  four.  Nonlinear  response  of  these  control  actuators  is  simulated  by 
saturation  and  rate  limiting  the  output  torques  as  shown  in  Figure  11.  Wheel  torque 
commands  are  generated  by  combining  feed  forward  and  error  based  control  laws.  The 
control  laws  implementation  is  delayed  at  initiation  to  allow  the  attitude  determination 
system  to  converge.  High  gain  control  laws  dominate  the  wheel  torque  response.  Figures 
43  and  44  show  the  torque  and  momentum  response  of  the  three  operating  reaction 
wheels. 
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Wheel  Torques 
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Figure  43.  Reaction  Wheel  Control  Torques 


Attitude  control  performance  is  dominated  by  attitude  determination  errors.  The 

control  system  is  designed  to  limit  steady  state  errors  in  the  attitude  quaternions  to  less 

than  2xl0"7  with  ideal  attitude  and  angular  rate  knowledge.  With  the  control  laws 

applied,  the  attitude  control  response  resembles  the  errors  of  the  attitude  determination 

system  in  steady  state.  Figure  45  shows  the  attitude  quaternion  error  response  for  the 

baseline  simulation  run.  Damped  corrections  back  to  steady  state  after  the  star  gap  are 

much  slower  than  those  realized  by  the  attitude  determination  system. 
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Figure  45.  Baseline  Control  Attitude  Quaternion  Error 


The  error  in  the  controlled  spacecraft  angular  rate  is  shown  in  Figure  46.  Since 
most  of  the  attitude  determination  noise  is  removed  by  the  controller,  the  uncorrected  rate 
error  bias  due  to  the  star  gap  is  evident.  Although  the  errors  remain  small  during  the  star 
gap,  the  bias  drives  the  attitude  quaternion  error  build  up. 


Figure  46.  Baseline  Control  Angular  Rate  Error 


B.  DYNAMIC  GYRO  VS  GYRO  PERFORMANCE 

Results  presented  in  this  section  compare  a  gyro  based  attitude  determination 
system  with  the  dynamic  gyro  based  system.  Three  orthogonally  mounted  gyroscopes 
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are  modeled  to  provide  simulated  spacecraft  body  angular  rate  measurements  to  the 
attitude  propagator  when  the  gyro  option  is  selected.  The  gyro  model  includes  static  bias, 
rate  noise,  and  rate  random  walk  supplied  by  integration  of  random  noise.  Gyroscope 
characteristics  used  for  the  simulation  results  are  listed  in  Table  3.  A  bias  error  tracking 
system  is  also  modeled  in  parallel  with  the  gyro  generated  rates  to  accept  bias  corrections 
provided  by  the  Kalman  filter.  Direct  comparison  of  the  dynamic  gyro  and  gyro  based 
attitude  determination  systems  is  accomplished  using  identical  input  parameters  to  obtain 
simulation  results.  The  SIMULINK  subsystem  for  gyro  rate  measurement  and  bias  error 
tracking  is  illustrated  in  Figure  47. 


Figure  47.  SIMULINK  Subsystem  Diagram:  Gyro  Measurement  and  Bias  Error 

Tracker. 

1.  Continuous  Star  Tracker  Coverage 

With  precise  attitude  updates  from  star  trackers  there  is  only  little  noticeable 
difference  between  the  dynamic  gyro  and  gyro  based  attitude  determination  systems. 
Figure  48  provides  a  side-by-side  comparison  of  simulation  results  obtained  with  a 
continuous  star  update  interval  of  2  seconds. 


76 


Figure  48.  Results:  Dynamic  Gyro  vs.  Gyro  with  Continuous  Star  Coverage 


Although  the  attitude  quaternion  errors  show  little  difference  between  the  two 
systems  the  rate  errors  from  the  gyros  are  noticeably  more  accurate  and  better  resemble 
white  noise.  Since  the  characteristics  of  the  rate  errors  are  different  between  the  two 
systems  the  plant  error  covariance  matricies  used  in  the  Kalman  filter  are  chosen 
differently.  Table  4  shows  the  plant  error  covariance  matricies  used  by  the  dynamic  gyro 
and  gyro  based  systems.  These  values  were  determined  through  tuning  over  several 
simulation  runs. 
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Plant  Error  Covariance  Matrix  (Q) 

Dynamic  Gyro 

Gyro 
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0.5 
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0 

0.03 

0 

_0 
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0 

0 

0.5_ 

_0 

0 

0 

0 

0 

0.03_ 

Table  4.  Plant  Error  Covariance  Matricies 


2.  Gapped  Star  Tracker  Coverage 

Figure  49  shows  the  comparison  of  the  dynamic  gyro  and  gyro  based  attitude 
determination  systems  when  a  200  second  gap  in  star  tracker  coverage  is  encountered 
during  maneuvering  operations.  Very  little  difference  is  evident  in  the  rate  error  plots 
because  the  bias  build  up  during  the  star  gap  is  so  small  compared  to  the  noise  in  the 
error.  The  effect  of  the  bias,  however,  shows  in  the  attitude  quaternion  error  plots.  At 
the  end  of  the  star  gap  the  error  in  the  dynamic  gyro  based  system  is  about  five  times  that 
of  the  gyro  based  system. 
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DG:  Quaternion  Errors 
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Figure  49.  Results:  Dynamic  Gyro  vs.  Gyro  with  Gapped  Star  Coverage 


C.  DYNAMIC  GYRO  PLANT  ERROR  ANALYSIS 


In  this  section  the  effects  of  major  error  sources  on  the  performance  of  the 
dynamic  gyro  are  analyzed.  Figure  50  shows  the  attitude  determination  plots  for  the 
baseline  simulation  run.  Dynamic  gyro  input  parameters  are  varied  and  results  are 
compared  with  these  plots. 
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Figure  50.  Baseline  Attitude  Determination  Performance  Results 


1.  Disturbance  Torque  Modeling 

The  accuracy  of  the  dynamic  gyro  depends  on  knowledge  of  external  disturbance 
moments  and  internal  spacecraft  momentum.  For  the  baseline  simulation  run  no  external 
disturbances  were  modeled  in  the  dynamic  gyro.  This  leads  to  the  relatively  large  drift  in 
attitude  quaternions  when  uncorrected  during  the  star  gaps.  If  the  spacecraft’s  attitude 
with  respect  to  the  orbital  reference  frame  can  be  determined  the  gravity  gradient 
disturbance  torque  can  be  modeled.  For  the  Bifocal  Relay  Mirror  satellite  this 
disturbance  has  the  greatest  effect.  Figure  51  illustrates  the  increase  in  performance  of 
the  dynamic  gyro  based  attitude  determination  system  when  the  gravity  gradient  moment 
is  modeled.  At  the  end  of  the  star  gap  the  attitude  errors  are  comparable  to  the  gyro 
based  system.  The  total  system  angular  momentum  error  is  much  smaller  during  the  star 
gap  and  improvement  can  even  be  seen  during  times  of  continuous  star  coverage. 
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Figure  51.  Gravity  Gradient  Disturbance  Modeled  in  the  Dynamic  Gyro 


2.  Rotation  Axis  Alignment  Error 

Alignment  errors  of  momentum  exchange  control  devices  and  slewing 
appendages  have  direct  effects  on  the  dynamic  gyro  momentum  error.  Alignment  errors 
of  fixed  reaction  wheels  are  easily  compensated  for  by  calibration  but  control  moment 
gyros  have  directionally  variant  momentum  vectors  with  respect  to  the  spacecraft  body. 
The  results  plots  shown  in  Figure  52  are  produced  when  a  periodic  alignment  error  on  the 
net  momentum  of  the  reaction  wheels  with  a  magnitude  of  approximately  0.5  degrees  is 
added  to  the  simulation.  A  significant  increase  in  attitude  error  is  developed  during  the 
star  gap.  The  dynamic  gyro  does  not  track  the  system  angular  momentum  as  well  even 
during  continuous  star  coverage. 
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Figure  52.  Reaction  Wheel  Alignment  Error  Effects  on  the  Dynamic  Gyro 
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3. 


Potentiometer  Quantization 


The  patterns  shown  in  angular  rate  error  plots  for  the  dynamic  gyro  do  not  appear 
as  white  noise  during  appendage  motion  because  of  quantization  effects  in  the  model  of 
the  potentiometer  that  measures  the  relative  orientation.  The  simulated  quantization  level 
is  0.01  degrees  in  the  baseline  simulation.  Since  the  effect  alternates  direction  of  error 
the  attitude  quaternion  errors  are  not  affected  significantly.  Figure  53  shows  simulation 
performance  results  when  the  quantization  level  is  decreased  an  order  of  magnitude.  The 
angular  rate  error  looks  more  like  white  noise  and  the  bias  build  up  during  the  star  gap 
can  be  seen. 
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Figure  53.  Reduced  Potentiometer  Quantization  Effect  on  Dynamic  Gyro 

Performance 


4.  Moment  of  Inertia  Update  F requency 

The  inertia  matrix  for  spacecraft  with  relatively  small  or  slowly  moving 
appendages  does  not  change  quickly.  In  these  cases  processing  power  may  be  saved  by 
performing  the  system  inertia  calculation  at  a  lower  bandwidth  than  the  dynamic  gyro.  In 
the  Bifocal  Relay  Mirror  satellite  the  secondary  body  is  large  so  even  small  slew  rates 
cause  significant  change  in  the  system  inertia  matrix.  Figure  54  shows  the  effects  of 
decreasing  the  inertia  calculation  frequency  from  the  20  Hz  rate  of  the  dynamic  gyro  to 
once  every  10  seconds.  The  quaternion  error  profile  is  significantly  altered  but  the 
magnitude  of  the  error  is  only  slightly  increased.  The  momentum  error  in  the  dynamic 
gyro  takes  longer  to  correct  after  the  star  gap  since  the  Kalman  filter  attempts  to  correct 
for  all  errors  as  if  they  were  due  to  spacecraft  momentum. 
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Figure  54.  10  Hz  Moment  of  Inertia  Calculation 

D.  STAR  TRACKER  MEASUREMENT  ERROR  ANALYSIS 

The  dynamic  gyro  based  attitude  determination  system  is  highly  reliant  on 
Kalman  filter  updates  based  on  star  tracker  measurements.  The  quality  of  the  Kalman 
filter  corrections  depend  on  the  accuracy  of  the  trackers,  the  discrete  measurement 
interval  and  number  of  trackers  providing  measurements. 

1.  Star  Tracker  Accuracy 

The  results  shown  in  Figure  55  are  from  simulations  using  random  selection  at  2 
second  intervals.  The  plot  on  the  left  is  generated  using  a  1  sigma  variance  of  lxl0"4in 
the  horizontal  and  vertical  measurements  of  the  trackers.  To  generate  the  plot  on  the 
right  the  variance  is  reduced  by  four  times.  The  direct  effect  on  attitude  determination 
performance  is  apparent. 
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Figure  55.  Star  Tracker  1  Sigma  Variance:  lxl 0"4  [left]  vs  2.5x1  O'5  [right] 

2.  Update  Interval 

Figure  56  shows  a  comparison  of  attitude  determination  performance  using  2  and 
6  second  star  measurement  intervals.  Shorter  update  intervals  result  in  quicker 
convergence  and  less  steady  state  error. 
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Figure  56.  Star  Tracker  Update  Interval:  2  Seconds  [left]  vs  6  Seconds  [right] 

3.  Star  Tracker  Selection 

Because  the  star  tracker  measurements  provide  useful  information  along  only  two 
axes  at  least  two  trackers  must  be  active  to  maintain  pointing  control.  Figure  57  shows 
the  comparison  of  an  attitude  determination  simulation  with  updates  spread  evenly 
between  three  trackers  and  one  that  uses  two  trackers  with  95%  of  the  updates  coming 
from  the  same  sensor.  There  is  no  apparent  degradation  when  the  updates  are  spread 


evenly  between  the  two  operational  trackers. 
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Figure  57.  Star  Tracker  Selection  Comparison:  Even  Distribution  [left]  vs  One 

Tracker  Favored  [right] 
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VIII.  CONCLUSIONS 


Dynamic  modeling  provides  an  imperfect  but  operative  means  of  estimating 
multi-body  spacecraft  angular  rates  when  mechanical  gyro  data  are  not  available.  An 
attitude  and  angular  rate  estimation  scheme  is  developed  in  this  thesis  that  integrates  the 
dynamic  gyro  concept  with  an  error  state  extended  Kalman  filter  estimator  that  provides 
precise  attitude  updates  from  star  tracker  measurements.  Results  indicate  that  the 
determination  system  provides  effective  estimates  for  performing  attitude  control. 

A.  SUMMARY 

The  attitude  determination  system  is  incorporated  into  a  multi-body  spacecraft 
attitude  simulation  for  evaluation  and  analysis.  Simulation  is  an  extremely  valuable 
analysis  tool  for  understanding  the  effects  of  error  sources  on  system  performance.  It  is 
also  a  suitable  mechanism  for  comparison  with  gyro  based  determination  systems  since  a 
mechanical  gyro  model  can  easily  be  inserted. 

The  effects  of  the  primary  error  sources  on  the  dynamic  gyro  performance  are 
investigated  through  simulation.  It  is  shown  that  the  corrections  provided  by  a  star 
tracker  based  Kalman  filter  make  the  system  robust  to  measurement  and  parameter 
knowledge  error  sources.  Significant  improvement  in  attitude  determination  performance 
is  realized  when  disturbance  torques  are  modeled.  The  other  primary  error  sources 
include  the  alignment  of  momentum  exchange  control  devices  and  relative  angle  and  rate 
knowledge  in  large  or  quickly  slewing  appendages.  Error  effects  are  amplified  during 
star  gaps  when  no  corrections  to  the  dynamic  model  are  available. 

The  software  implemented  dynamic  gyro  essentially  emulates  the  functions  of  a 
set  hardware  gyroscopes.  In  a  spacecraft  where  the  mechanical  gyros  have  failed  or 
become  too  erratic  to  be  corrected  by  the  Kalman  filter  the  dynamic  gyro  may  be  a  viable 
replacement.  Operated  in  tandem  with  mechanical  gyros,  either  system  can  provide 
redundant  inertial  angular  velocity  for  improved  attitude  control. 
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This  attitude  determination  concept  is  ideally  suited  to  a  spacecraft  designed 
specifically  for  its  implementation  with  precise  internal  sensors  and  mechanisms  to 
monitor  spacecraft  parameters  and  integrated  external  torque  estimation  modeling.  In 
these  systems  the  dynamic  gyro  can  increase  the  lifetime  and  reliability  of  the  spacecraft 
while  reducing  the  power  requirements. 

B.  RECOMMENDATIONS 

A  strong  recommendation  for  future  work  is  simulation  model  tailoring  to  real 
spacecraft  hardware.  The  multi-body  dynamics  can  be  expanded  to  model  an  actual 
satellite’s  mass  and  inertia  characteristics.  Additionally,  sensor  parameters  can  be 
modeled  to  match  existing  sensor  error  specs.  With  these  modifications  the  simulation 
model  can  be  used  conduct  analysis  and  predict  performance  when  the  dynamic  gyro 
software  is  implemented  on  board  the  spacecraft.  Future  work  may  include  using  the 
simulation  to  compare  the  predicted  dynamic  gyro  based  attitude  determination 
performance  to  actual  gyro  based  performance  recorded  with  telemetry  data. 

There  are  several  improvements  that  can  be  made  to  the  attitude  determination 
simulation  model.  It  could  easily  be  modified  to  provide  the  capability  to  mix  gyroscope 
and  dynamic  gyro  data  for  redundant  information.  A  parity  matrix  developed  from  the 
pseudo-inverse  concept  can  be  generated  to  account  for  system  observability  in  the  over 
determined  case.  This  would  allow  performance  analysis  with  selected  gyro  failures  and 
further  aid  in  the  evaluation  of  its  utility  as  a  back-up  attitude  determination  scheme. 
Another  simulation  improvement  would  be  the  incorporation  of  higher  order  dynamic 
effects  into  the  model  including  center  of  mass  offsets  and  flexibility  modes  for  the 
appendage  couplings. 

Finally,  improvements  to  the  attitude  determination  scheme  can  be  developed 
including  a  torque  error  estimator  as  suggested  by  The  Aerospace  Corporation  [Ref.  8]. 
The  better  the  dynamics  are  modeled  in  the  processing  software  the  more  effective  the 
dynamic  gyro  becomes  as  a  replacement  for  the  hardware  gyroscope. 
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APPENDIX  A:  KALMAN  FILTER  BACKGROUND 


“The  Kalman  filter  combines  all  available  measurement  data,  plus  prior 
knowledge  about  the  system  and  measuring  devices,  to  produce  an  estimate  of  the  desired 
variables  in  such  a  manner  that  the  error  is  minimized  statistically”  [Ref.  2].  An 
understanding  of  the  Kalman  filter  requires  some  background  in  the  theory  of  probability 
of  random  variables  and  processes  [Ref.  1], 


A.  PROBABILITY 


The  probability  of  an  event,  e,  represents  a  possible  outcome  of  a  random 
experiment  and  is  written  Pr(e).  A  random  variable,  X,  can  be  thought  of  as  a  function  of 
the  outcomes  of  some  random  experiment.  The  manner  of  specifying  the  probability  with 
which  different  values  are  taken  by  a  random  variable  is  the  probability  distribution 
function,  F(x) 


F(x)=Pr(X  <  x)  (A.l) 

This  is  a  function  over  the  range  of  possible  values  that  shows  the  probability  with  which 
the  random  variable  takes  on  a  value  at  or  less  than  the  value  of  the  range.  Its  derivative 
is  the  probability  density  function,  f(x), 


f(x)= 


dF(x) 

dx 


(A.2) 


This  function  identifies  the  likelihood  of  a  random  variable  assuming  a  particular  value  in 
its  range  of  possible  values.  Over  the  range  of  all  possible  values  a  characteristic  of  any 
probability  distribution  or  density  function  is 


F(oo)=Jf(u)du=l  (A.3) 

A  joint  probability  density  function  can  be  defined  for  multiple  variables.  For  two 
random  variables  the  joint  density  is  given  by 


f2(x>y) 


52f2(xa) 

dxdy 


(A.4) 
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B.  EXPECTATION 


The  expectation  of  a  random  variable  is  defined  as  the  sum  of  all  values  the 
random  variable  may  take,  each  weighted  by  the  probability  with  which  the  value  is 
taken. 


E[X]=|  xf(x)dx  (A.5) 

-oo 

which  is  also  called  the  mean  value  or  first  moment  of  X.  This  is  a  precisely  defined 
number  toward  which  the  average  of  a  number  of  observations  of  the  random  variable  X 
tends.  Since  a  function  of  a  random  variable  is  itself  a  random  variable,  the  expectation 
of  a  function  of  X,  E[g(X)],  can  be  expressed  as 

oo 

E[g(X)]= j  g(x)f(x)dx  (A.6) 

An  important  statistical  parameter  descriptive  of  X  is  its  mean  squared  value  defined  by 


E[X2]  =  J  x2f(x)dx 


(A.7) 


which  is  also  called  the  second  moment  of  X.  The  root- mean- squared  (rms)  value  of  X  is 
-\/e[X2]  .  The  variance  of  a  random  variable  is  defined  as  the  mean  squared  deviation  of 
the  variable  from  its  mean  denoted  by  6 2 . 


o2  =  J  (x-E[X])2f(x)dx  =  E[X2  ]-E[X]2  (A.8) 

For  zero  mean  random  variables  the  variance  is  simply  E[X2] .  The  square  root  of  the 
variance,  6  ,  is  called  the  standard  deviation  of  the  random  variable. 

A  very  important  concept  is  that  of  statistical  correlation  between  random 
variables.  A  partial  indication  of  the  degree  to  which  one  variable  is  related  to  another  is 
given  by  the  cross  covariance,  which  is  the  expectation  of  the  product  of  the  deviations  of 
the  two  variables  from  their  means, 
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E[(X-E[X])(Y-E[Y])]=  1 1  (x-E[X])(y-E[Y])f2(x,y)dydx=E[XY]-E[X]E[Y] 

(A.9) 

For  a  vector  of  random  variables  a  symmetric  covariance  matrix  can  be  defined 
where  the  diagonal  elements  are  the  individual  variances  of  the  vector  components  and 
the  off  diagonals  are  given  by  the  cross  covariances  between  the  corresponding  vector 
components.  The  cross  covariance,  normalized  by  the  standard  deviations  of  X  and  Y,  is 
called  the  correlation  coefficient. 

ft=EMX]E[Y]  (A10) 

6x6y 

This  is  a  measure  of  the  degree  of  linear  dependence  between  variables  X  and  Y.  If  they 
are  completely  independent  n  is  zero. 

C.  LEAST-SQUARES  ESTIMATION  CONCEPT 

The  optimality  condition  used  in  the  Kalman  filter  is  the  minimization  of 
weighted  least- squares  error.  The  least-squares  minimization  problem  involves  a  set  of 
measurements,  y ,  which  are  linearly  related  to  the  vector  of  state  variables,  x ,  by  the 
expression 

y=Hx+v  (A.ll) 

where  v  is  an  unknown  vector  of  actual  measurement  errors  with  zero  mean.  The  error 
we  seek  to  minimize  is  the  measurement  residual,  z ,  given  by 

z=y-Hx  (A.  12) 

where  x  is  the  estimate  of  the  actual  state  vector.  Since  the  sum  of  the  squares  of  a 
vector  are  generated  by  the  inner  product,  the  cost  function,  J,  is 

J=(y-Hk)T(y-Hk)  (A.  13) 

Minimization  is  obtained  by  differentiating  and  setting  the  result  to  zero 
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(A.  14) 


x=(hth)  '  HTy  (A.  16) 

A  weighted  least- squares  solution  can  be  used  when  not  all  components  of  the 
measurement  residual  are  treated  equally 

x=(htWH)  '  HTWy  (A.  17) 

where  W  is  a  positive  semidefinite  matrix  relating  scale  factors  between  components. 
Expanding  this  solution  to  multiple  estimates  over  time  assumes  all  measurements  are 
used  together  in  a  batch  processing  scheme.  Every  time  a  new  measurement  becomes 
available  it  is  appended  to  the  measurement  vector  and  the  estimated  state  vector  includes 
estimates  corresponding  to  each  of  the  accumulated  measurements.  Hie  Kalman  filter  is 
based  on  recursive  processing  where  each  measurement  is  used  sequentially  to  generate 
an  optimal  estimate  of  the  current  state  without  recomputing  estimates  of  all  previous 
states.  For  this  technique  all  previous  information  is  embodied  in  the  prior  estimate  and 
state  covariance  matrix. 


D.  STATE  ESTIMATE  AND  COVARIANCE 

The  Kalman  filter  estimation  algorithm  maintains  the  first  two  statistical  moments 

of  the  estimated  state.  The  estimated  state,  x ,  is  a  vector  of  random  variables  whose 
mean  (first  moment)  is  the  actual  state  vector,  x .  The  error  in  the  state  estimate  is 
defined  as 

xDx-x  (A.18) 

and  is  thus  assumed  to  be  zero  mean.  The  covariance  of  the  state  estimate  error, 
designated  P,  is  given  by 
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P=E[Mt]  (A.  19) 

This  covariance  matrix  associated  with  the  state  estimate  error  (second  moment)  provides 

a  statistical  measure  of  the  uncertainty  in  x  and  the  correlation  between  the  errors  of  its 
components.  For  a  state  vector  of  two  variables  the  state  estimate  error  is  given  by 


and  has  the  covariance  matrix 


(A.20) 


IT  *1  VG 

L 

"  E[x2] 

E[X[  x2  ] 

lUi*2  *2  _ 

J 

EjXjXj 

E[x22]  _ 

(A.21) 


The  diagonal  elements  of  the  state  covariance  matrix  are  the  mean  square  errors  in 
the  knowledge  of  the  state  variables  while  the  off  diagonal  elements  are  indicators  of 
cross  correlation  between  elements  of  the  estimated  error. 


E.  STOCHASTIC  LINEAR  DYNAMIC  MODEL 

The  Kalman  filter  requires  representation  of  system  dynamics  in  a  linearized 
state-space  form,  a  linear  measurement  model,  and  assumed  characteristics  of  process 
and  measurement  noise.  For  a  continuous  linear  system  the  general  state-space  model 
and  measurement  equations  are  given  by 

x(t)=F(t)x+G(t)w(t) 

y(t)=H(t)x(t)+v(t) 

where  w(t)  and  v(t)  are  random  vectors  representing  the  unmodeled  disturbance  inputs 
and  measurement  errors.  These  vectors  are  treated  as  unbiased  (zero  mean)  white  noise. 
Note  that  if  one  of  these  vectors  is  known  to  have  a  nonzero  bias  the  mean  can  be 
augmented  onto  the  state  vector  creating  a  new  state  space  model  with  unbiased  random 
error.  The  covariance  matricies  associated  with  the  process  disturbance,  w(t) ,  and  the 
measurement  noise,  v(t) ,  are 
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Q  =  E[v(t)v(t)T] 
R  =  E[w(t)w(t)T] 


(A.23) 


The  equivalent  discrete  representation  is 


(A.24) 


x9+i^ ¥  k  k 
yk=Hkxk+Vk 

Here  the  subscript  k  represents  values  at  time  k  and  subscript  k+1  represents  values  at  the 
next  discrete  time  step  tk+1.  Note  that  for  time  invariant  systems  the  discrete  state 
transition  matrix,  ok ,  is  related  to  the  continuous  formulation  by 


ok=ek+lk  (A.25) 

which  depends  only  on  the  systems  dynamics  matrix,  F,  and  the  discrete  time  interval.  If 
the  discrete  interval  is  kept  short  enough  the  relation  holds  as  an  approximation.  The 
state  transition  matrix  allows  the  calculation  of  state  vector  at  the  next  discrete  time  step 
in  the  absence  of  forcing  functions.  It  obeys  the  differential  equation 


do, 

dt 


=  F(t)ok 


(A.26) 


F.  PROPAGATION  OF  ERRORS 

The  recursive  Kalman  filter  requires  the  propagation  estimate  and  error 
covariance  based  on  system  dynamics.  In  the  discrete  implementation  the  error  in  the 
current  estimate  given  by 

xk=xk-xk  (A.27) 

has  the  covariance  matrix  representing  the  uncertainty  in  the  estimate 

Pk  =  E[xkxk]  (A.28) 

The  best  estimate  of  the  future  state,  xk+1 ,  is  given  by 

Xk+1=  k^k  (A.29) 

The  error  in  the  new  estimate 

xk+l=okxk-\wk  (A.30) 


has  the  expected  value 
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E[xk+1]  =  ok  E[xk  ]-  4  E[wk  ]=  0  (A.31) 

since  xk  and  wk  are  assumed  to  be  unbiased.  Thus  the  predicted  estimate  remains 
unbiased.  Note  that  if  a  deterministic  input  is  added  to  the  dynamic  system  model  then 
the  identical  quantity  is  added  to  both  the  actual  and  estimated  state  leaving  the 
estimation  error  unchanged.  It  can  be  shown  that  the  associated  state  error  covariance 
matrix  is  given  by 

Pk+i  =  E[xk+1xk+1l=  okPkok  +  4  Q  4  (A.32) 

where  Q;  is  the  covariance  of  the  random  system  disturbance.  It  is  evident  from  this 
equation  that  the  size  of  disturbance  directly  impacts  the  uncertainty  in  the  state  estimate. 
The  system  dynamic  stability  reflected  in  the  state  transition  matrix  also  effects  the 
uncertainty.  The  covariance  of  neutrally  stable  or  unstable  systems  will  grow  unbounded 
over  time  in  the  absence  of  measurement  corrections.  The  propagation  of  the  state 
estimate  and  its  error  covariance  matrix  is  the  prediction  step  of  the  Kalman  filter 
algorithm. 

G.  MEASUREMENT  UPDATES 

The  correction  step  of  the  Kalman  filter  algorithm  incorporates  measurement 
updates  to  refine  the  state  estimate  and  error  covariance.  This  step  is  executed  only  when 
a  measurement  becomes  available.  When  a  measurement  is  taken  we  use  the  superscripts 
and  ‘+’  to  denote  values  at  a  particular  time  before  and  after  incorporation  of  the 
measurement  correction.  Given  a  prior  estimate  of  the  system  state  at  time  k  we  seek  to 
update  our  estimate  based  on  the  measurement  yk  in  a  linear,  recursive  form 

xk=Kkik  +  Kkyk  (A.33) 

where  Kk  and  Kk  are  matricies,  as  yet  unspecified,  that  determine  the  relative  weighting 
of  the  prior  estimate  and  current  measurement.  The  error  in  this  estimate  can  be  shown  to 
be 

K  =  (<  +KtHt  -l)x„  +  Kk*k  +  Kkvk  (A.34) 

Since  xk  and  vk  are  unbiased 
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E[xk]=0 
E[vk]  =  0 

this  estimate  can  only  remain  unbiased  for  all  given  states  if 
Kk=I-KkHk 

With  this  requirement  the  estimator  takes  the  form 

x^=xk  +  Kk(yk-Hkkk) 
and  has  corresponding  error 

=  (I~KkHk)kk  +  Kkvk 


(A.35) 


(A.36) 


(A.37) 


(A.38) 


The  state  error  covariance  must  also  be  updated. 

Pk  =  E[xkxkT] 

=  E{(I-  KkHk)Xk[XkT  (I  -  KkHk)T  +  vkKk  ]  +  Kkvk[kkT  (I  -  KkHk)T  +  vkKk  ]} 


(A.39) 

Using  the  definitions 

pk  =  E[kkrkTJ  (A.40) 

and 

Rk  =  E[vkvk  ]  (A.41) 

and  noting  that  the  measurement  errors  are  uncorrelated 

E[kkvk  ]  =  E[vkkkT  ]  =  0  (A.42) 

The  covariance  can  be  simplified  to 

Pk+  =  (I-  KkHk)Pk  (I-  KkHk)T  +  KkRkKk  (A.43) 

The  criterion  for  choosing  the  optimal  Kk  is  to  minimize  a  weighted  scalar  sum  of 
the  diagonal  elements  of  the  error  covariance  matrix,  Pk  .  Thus  the  cost  function  is 

lk  =  E[^Srk]  (A.44) 

where  S  is  any  positive  semidefinite  matrix.  Choosing  S  =  I, 
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Jk=trace[Pk+]  (A.45) 

which  is  equivalent  to  minimizing  the  length  of  the  estimation  error  vector.  To  find  the 
value  of  Kk  which  yields  a  minimum,  it  is  necessary  to  take  the  derivative  of  4  with 
respect  to  Kk  and  set  it  equal  to  zero.  The  result  is 

-2(1- KkHk)PkHk  +2KkRk  =0  (A.46) 

Solving  for  Kk, 

Kk=PkH^[HkPk^r1  (A.47) 

which  is  referred  to  as  the  Kalman  gain  matrix.  Since  the  Hessian  of  Jk  is  positive 
semidefinite  so  Kk  does  indeed  produce  a  minimum.  Substituting  into  the  equation  for 
the  updated  error  covariance  gives 

Pk*  =  Pk-PkH;[H^lH;+Rkr'H11Pk  =  (I-KlH1)Pk  (A.48) 
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APPENDIX  B:  MATLAB  CODE 


A.  INPUT  PARAMETERS  AND  SIMULATION  CALL 

%Bifocal  Relay  Mirror  Spacecraft  Attitude  Control  System 
%Input  parameters  and  control  options  for  attitude  simulation  model 

clear 

tic 

stoptime=500;  %Simulation  stoptime 
dt=.05;  %fixed  step  size  and  computer  time  step 

ml=2267.6059;  %Mass  of  X-mit  telescope  (kg) 
m2=972.3628;  %Mass  of  RCV  telescope  (kg) 

%Orbital  parameters  (circular) 

Re=6378.1363;  %Earth  radius  (km) 

mu=3.986004415e5;  %gravitational  constant  (kmA3/sA2) 

h=715;  %orbit  altitude  (km) 

r=Re+h;  %orbit  radius  (km) 

we=7.2921158533e-5;  %Earth  rotation  rate  (r/s) 

wO=sqrt(mu/rA3);  %orbital  rate  based  on  circular  orbit/altitude  (r/s) 

nu0=0;  %initial  true  anomaly  (r) 

u0=0;  %initial  angle  of  the  magnetic  dipole  normal  to  the 
%line  of  ascending  nodes  (r) 
el=ll*pi/180;  %magnetic  dipole  tilt  angle 
inc=40*pi/180;  %inclination 
Km=7.943el5;  %magnetic  constant  (NmA2/a) 

Bm=Km/(r*1000)A3;  %magnetic  field  strength  (N/am) 

%Disturbance  torques 

Mds=le-4*[1;1;1];  %  secular  disturbance  torques  (Nm) 
Mdp=4e-4*[1,1,1];  %periodic  disturbances  torque  magnitudes  (Nm) 
pdper=[400, 500,600];  %periodic  disturbance  torque  periods  (s) 
pdfreq=2*pi./pdper;  %periodic  disturbance  torque  frequencies  (r/s) 
pdphase=[pi/4, pi/2, pi];  %periodic  disturbance  torque  phase  (r) 

%Magnetic  dumping  control 

magoo=0;  %magnetic  dumping  on/off  control  (0=off,  l=on) 
magsat=180;  %maximum  torque  rod  output 
magk=5e5;  %magnetic  torquer  gain 

%Appendage  (RCV  scope)  motion  relative  to  body  x-axis  (sinusoidal) 
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bper=300;  %period  of  oscillation  (s) 
bstart=.25*bper;  %start  motion  (s) 
bmotion=l*bper;  %length  of  motion  (s) 
bstop=bstart+bmotion;  %stop  motion  (s) 
bfreq=2*pi/bper;  %frequency  of  motion  (rad/s) 
bamp=15*pi/180;  %amplitude/half  maximum  angular  offset  (rad) 
bi=0;  %initial  angular  offset  (rad) 
bdi=0;  %initial  relative  angular  rate  (rad/s) 

bacc=bamp*bfreqA2;  %maximum  realtive  angular  acceleration  (rad/sA2) 
bnv=le-5A2;  %relative  angular  noise  variance  from  potentiometer  (rad) 
bns=0;  %initial  seed  of  potentiometer  noise 
bq=.01*pi/180;  %quantization  of  potentiometer  readings  (rad) 

%Command  profile  inertial  to  body  accelerations  (sinusoidal) 
wdcfreq=2*pi*[l/800;l/500;l/600];  %frequency  of  acceleration  profile  variation 
wdc=[.05*wdcfreq(l)A2;.5*wdcfreq(2)A2;-.l*wdcfreq(3)A2]; 

%amplitude  of  acceleration  profile 
wdcphase=[pi/8;-pi/12;0];  %phase  of  acceleration  profile 
cpoo=l;  %tum  commanded  profile  on/off  (0=off,  l=on) 

Il=[2997.28025, -3.9331, 118.2824; 

-3.9331,3164.18285,1.1230; 

118.2824,1.1230,881.82105]; 

%moment  of  inertia  matrix  of  primary  body  (XMIT  scope)  about  its  cm 
I2=[1721.07340,-.01 16,-7.8530; 

-.0116,1559.85414,-12.5463; 

-7.8530,-12.5463,182.89962]; 

%moment  of  inertia  matrix  of  appendage  (RCV  scope)  about  its  cm 

dxl=.558354158;  %Distance  of  X-mit  telescope  cm  from  S/C  cm  in  x  dir  (m) 
dyl=3.91788e-4;  %Distance  of  X-mit  telescope  cm  from  S/C  cm  in  y  dir  (m) 
dzl=-.  15226902;  %Distance  of  X-mit  telescope  cm  from  S/C  cm  in  z  dir  (m) 
dx2=- 1.302 1139 18;  %Distance  of  RCV  telescope  cm  from  S/C  cm  in  x  dir  (m) 
dy2=-9. 13673e-4;  %Distance  of  RCV  telescope  cm  from  S/C  cm  in  y  dir  (m) 
dz2=. 355 100092;  %Distance  of  RCV  telescope  cm  from  S/C  cm  in  z  dir  (m) 

IC 1  =1 1 +m  1  *  [dy  1  A2+dz  1 A2  ,dx  1  *dy  1  ,dx  1  *dz  1 ; 
dx  1  *dy  1  ,dx  1  A2+dz  1  A2,dy  1  *dz  1 ; 
dx  1  *dz  1  ,dy  1  *dz  1  ,dx  1  A2+dy  1 A2] ; 

%moment  of  inertia  matrix  of  primary  body  (XMIT  scope)  about  S/C  cm 
roti=[l,0,0;0,cos(bi),sin(bi);0,-sin(bi),cos(bi)]; 

%rotation  matrix  corresponding  to  initial  relative  angular  position 
%of  appendage  (RCV  scope) 

%x-axis  rotation  of  magnitude  bi 

Ic2i=roti'*I2*roti;  %initial  rotated  moment  of  inertia  matrix  of  appendage 
%(RCV  scope) 
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I2m=m2*[dy2A2+dz2A2,dx2*dy2,dx2*dz2; 

dx2*dy2,dx2A2+dz2A2,dy2*dz2; 

dx2*dz2,dy2*dz2,dx2A2+dy2A2]; 

%moment  of  inertia  of  appendage  (RCV  scope)  about  S/C  C.M. 
%due  to  mass  offset 

IC2i=roti'*I2*roti+I2m;  %initial  moment  of  inertia  matrix  of  appendage 
%(RCV  scope)  about  S/C  cm 
Ii=ICl+IC2i;  %initial  S/C  moment  of  inertia  about  S/C  cm 

%Wheel  control  law  gains 

%u=-Kqqe-Kwwe=-Kq(qc*q)-Kw(w-wc) 

cldel=30;  %delay  control  laws  for  attitude  determination 

Kql=3000; 

Kq2=7000; 

Kq3=4500; 

Kwl=1000; 

Kw2=2000; 

Kw3=1000; 

%Initial  body  with  respect  to  inertial  quaternions 
qli=0;q2i=0;q3i=0; 

qi=[qli,q2i,q3i,sqrt(l-qliA2-q2iA2-q3iA2)]; 

%Initial  body  with  respect  to  orbit  quaternions 
qloi=0;q2oi=0;q3oi=0; 

qoi=[qloi,q2oi,q3oi,sqrt(l-qloiA2-q2oiA2-q3oiA2)]; 

%Initial  body  with  respect  to  inertial  direction  cosine  matrix 

dcmi(  1 ,  l)=qi(  1 )  A2-qi(2)A2  -qi(3)  A2+qi(4)A2; 

dcmi(  1 ,2)=2*(qi(  1  )*qi(2)+qi(3)*qi(4)) ; 

dcmi(  1 ,3)=2*(qi(  1  )*qi(3)-qi(2)*qi(4)) ; 

dcmi(2,l)=2*(qi(l)*qi(2)-qi(3)*qi(4)); 

dcmi(2,2)=-qi(l)A2+qi(2)A2-qi(3)A2+qi(4)A2; 

dcmi(2,3)=2*(qi(2)*qi(3)+qi(l)*qi(4)); 

dcmi(3, 1  )=2*(qi(  l)*qi(3)+qi(2)*qi(4)); 

dcmi(3,2)=2*(qi(2)*qi(3)-qi(l)*qi(4)); 

dcmi(3,3)=-qi(l)A2-qi(2)A2+qi(3)A2+qi(4)A2; 

Ci=[dcmi(  1,1)  ;dcmi(  1,2)  ;dcmi(  1,3); 
dcmi(2,l);dcmi(2,2);dcmi(2,3); 
dcmi(3 , 1 )  ;dcmi(3 ,2)  ;dcmi(3 ,3 )] ; 

%Direction  cosine  component  vector 
Ai=  [dcmi(  1,1)  ,dcmi(  1 ,2)  ,dcmi(  1,3); 
dcmi(2,l),dcmi(2,2),dcmi(2,3); 
dcmi(3,  l),dcmi(3,2),dcmi(3,3)] ; 

%Direction  cosine/attitude  matrix 
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^Initialize  S/C  momentum  with  wi  at  orbital  rate 

wi=[0;-w0;0];  %Initial  angular  rate  of  S/C  primary  body  due  to  orbital  rate 
Hi=Ii*wi;  %initial  angular  momentum  of  S/C  due  to  orbital  rate 

%Gyro  characteristics 

grsb=le-4*[-l, 1.5,1];  %arbitrary  gyro  static  bias 
grv=le-9*[l,l,l];  %variance  of  gyro  rate  noise 
gm=[0,l,2];  %initial  seed  of  gyro  rate  noise 
gasb=[0,0,0];  %gyro  acceleration  static  bias  (zeros) 
gav=le-12*[l,l,l];  %variance  of  gyro  acceleration  noise 
gan=[3,4,5];  %initial  seed  of  gyro  acceleration  noise 
grrwi=[0,0,0];  %initial  rate  random  walk  (zeros) 
wbi=[0;0;0];  %initial  gyro  bias  correction 
gx=0*pi/180;  %body  to  gyro(IRU)  x-axis  rotation 
gy=0*pi/180;  %body  to  gyro(IRU)  y-axis  rotation 
G=[cos(gy),sin(gy)*sin(gx),-sin(gy)*cos(gx); 

O,cos(gx),sin(gx); 

sin(gy),-cos(gy)*sin(gx),cos(gy)*cos(gx)]; 

%Gyro  alignment  matrix  -  body  to  IRU 

%Dynamic  Gyro  input  options 

dg=l;  %choose  Dynamic  Gyro  or  regular  gyros  (1=DG,  0=gyros) 
dggg=0;  %calculate  gravity  gradient  torque  for  DG  (l=yes(requires  or=l),  0=no) 
dgsd=0;  %secular  disturbance  torque  known  for  DG  (l=yes,  0=no) 
dgpd=0;  %periodic  disturbance  torque  known  for  DG  (l=yes,  0=no) 
moiu=l;  %steps  per  periodic  moment  of  inertia  calculation  (l=every  step) 
or=0;  %perform  orbital  reference  calculations 

%  Attitude  determination  initialization  with  error 
qeli=qli+.008;qe2i=q2i+.012;qe3i=q3i-.008; 
qei=[qeli,qe2i,qe3i,sqrt(l-qeliA2-qe2iA2-qe3iA2)]; 
dcmei(  1 , 1  )=qei(  1  )A2  -qei(2)  A2  -qei(3)  A2+qei(4)  A2 ; 
dcmei(  1 ,2)=2*(qei(  1  )*qei(2)+qei(3)*qei(4)); 
dcmei(  1 ,3)=2*(qei(  l)*qei(3)-qei(2)*qei(4)); 
dcmei(2, 1  )=2*(qei(  1  )*qei(2)-qei(3)*qei(4)); 
dcmei(2,2)=-qei(l)A2+qei(2)A2-qei(3)A2+qei(4)A2; 
dcmei(2,3)=2*(qei(2)*qei(3)+qei(l)*qei(4)); 
dcmei(3,l)=2*(qei(l)*qei(3)+qei(2)*qei(4)); 
dcmei(3,2)=2*(qei(2)*qei(3)-qei(l)*qei(4)); 
dcmei(3,3)=-qei(l)A2-qei(2)A2+qei(3)A2+qei(4)A2; 

Cei=[dcmei(  1 , 1 )  ;dcmei(  1,2); 
dcmei(l,3);dcmei(2,l);dcmei(2,2);dcmei(2,3); 
dcmei(3,l);dcmei(3,2);dcmei(3,3)]; 

Aei=[dcmei(  1 , 1  ),dcmei(  1 ,2),dcmei(  1,3); 
dcmei(2,l),dcmei(2,2),dcmei(2,3); 
dcmei(3,  l),dcmei(3,2),dcmei(3,3)] ; 
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qeloi=qloi+.02;qe2oi=q2oi+.03;qe3oi=q3oi-.02; 

qeoi=[qeloi,qe2oi,qe3oi,sqrt(l-qeloiA2-qe2oiA2-qe3oiA2)]; 

wei=wi+le-3*[-l;l;2]; 

Hei=Ii*wei; 

%Initialize  Kalman  filter 

xi=0*[ones(3,l)*eps*lel0;ones(3,l)*eps*le8];  %initial  state  vector 
Pi=[100*eye(3),zeros(3);zeros(3),eye(3)]*500;  %initial  covariance  matrix 
if  dg==l 
Qn=.5; 

else  %plant  noise  covariance  constant  DG  vs  gyro 

Qn=.03; 
end 

Q=[100*eye(3),zeros(3);zeros(3),eye(3)]*Qn;  %plant  noise  covariance 
%  small  Q  =  good  plant,  biq  Q  =  bad  plant 
R=10000*eye(2);  %measurement  noise  covariance  (10000) 

E=[1,0,0;0,1,0];  %choose  horizontal  and  vertical  components  only 

%Star  tracker  parameters 
corr=2*l/dt;  %computer  steps  between  stars 
gapstart=100;  %start  gap  in  star  tracker  data  (s) 
stargap=200;  %length  of  star  gap  (s) 
gapstop=gapstart+stargap;  %end  gap  in  star  tracker  data  (s) 
tlx=135*pi/180;  %body  to  star  tracker  one  x-axis  rotation 
tly=30*pi/180;  %body  to  star  tracker  one  y-axis  rotation 
Tl=[cos(tly),sin(tly)*sin(tlx),-sin(tly)*cos(tlx); 

0,cos(tlx),sin(tlx); 

sin(tly),-cos(tly)*sin(tlx),cos(tly)*cos(tlx)]; 

%body  to  star  tracker  one  rotation  matrix 
t2x=135*pi/180;  %body  to  star  tracker  two  x-axis  rotation 
t2y=-30*pi/180;  %body  to  star  tracker  two  y-axis  rotation 
T2=[cos(t2y),sin(t2y)*sin(t2x),-sin(t2y)*cos(t2x); 

0,cos(t2x),sin(t2x) ; 

sin(t2y),-cos(t2y)*sin(t2x),cos(t2y)*cos(t2x)]; 

%body  to  star  tracker  two  rotation  matrix 
t3x=180*pi/180;  %body  to  star  tracker  two  x-axis  rotation 
t3y=0*pi/180;  %body  to  star  tracker  two  y-axis  rotation 
T3=[cos(t3y),sin(t3y)*sin(t3x),-sin(t3y)*cos(t3x); 

0,cos(t3x),sin(t3x); 

sin(t3y),-cos(t3y)*sin(t3x),cos(t3y)*cos(t3x)]; 

%body  to  star  tracker  two  rotation  matrix 
mixT3=l;  %choose  whether  or  not  to  mix  in  star  tracker  3  data  (l=yes,  0=no) 
T3th=.3;  %choose  threshold  for  random  mix  of  star  tracker  3  data 
%(0=50%,  higher=less,  lower=more) 

T12b=0;  %choose  to  bias  random  selection  between  tracker  1  and  2 
%(0=50/50,  lower=favorl,  higher=favor2) 
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stnb=[0,0];  %unknown  star  tracker  bias  error 
stnv=le-4A2*[l,l];  %variance  of  star  tracker  noise 
stns=[0,l];  %initial  seed  of  star  tracker  noise 

%3  of  4  reaction  wheels  in  pyramid  constellation 
rwsat=l*[l,l,l];  %constellation  wheel  torque  saturation  limits  (Nm) 
rwrl=  10* [1,1,1];  %constellation  wheel  torque  rate  limiter  (Nm/s) 
hwi=-I2*[bdi;0;0];  %initialize  wheels  to  cancel  appendage  momentum 
rwa=45*pi/180;  %angle  of  reaction  wheels  to  x-y  plane 
RW=inv([-cos(rwa)*cos(pi/4),-cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4); 
-cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4); 
-sin(rwa),-sin(rwa),-sin(rwa)]);  %body  to  wheel  transform  matrix 
hwsi=RW*hwi;  %distribute  initial  wheel  momentum 
rwnm=[0;0;0];  %mean  of  wheel  noise 
rwnv=le-2A2*[l;l;l];  %variance  of  wheel  noise 
rwns=[7;8;9];  %initial  seed  of  wheel  noise 

rwev=le-8*[  1,1,1];  %variance  of  wheel  alignment  walk  error  noise 
rwes=[10, 11,12];  %initial  seed  of  gyro  acceleration  noise 
rwam=l*5e-3*[l,2,-l];  %reaction  wheel  constellation  alignment  error  magnitudes 
rwaf=2*pi*[  100,50,20];  %reaction  wheel  constellation  alignment  error  frequencies 
rwap=  [pi/4, pi/3, pi/2];  %reaction  wheel  constellation  alignment  error  phase 

sim('acs_sim509')  %Call  to  SIMULINK  simulation 

%Call  ploting  files  for  simulation  results  analysis 
profileplots  %command  profile 
ADplots509  %attitude  determination  performance 
ACplots509  %attitude  control  performance 

toe 


B.  SIMULATION  RESULTS  PLOTS 


1.  Commanded  Profile  Plots 

%Laser  Realy  Spacecraft  command  profile  plots 

figure(l) 

elf 

plot(tout,qc) 

title('Commanded  Quaternions') 
Iegend('qr,'q2','q3','q4') 
xlabel('Time  (sec)') 
grid  on 
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figure(2) 

elf 

plot(tout,wc) 

title('Commanded  Angular  Rates') 

Iegend('wl7w27w3') 

xlabel('Time  (sec)') 

ylabel('rad/s') 

grid  on 

figure(3) 

elf 

plot(tout,b*  180/pi) 
xlabel('Time  (sec)') 
ylabel('Deg') 

title('Relative  Angle  to  Recieve  Telescope  (alpha)') 

grid  on 

figure(4) 

elf 

Hm=sqrt(H( : ,  1 ) .  A2+H( :  ,2) .  A2+H( : ,  3) .  A2) ; 
plot(tout,Hm) 

title('System  Angular  Momentum') 
xlabel('Time  (sec)') 
ylabel('Nms') 
grid  on 


2.  Attitude  Determination  Performance  Results 

%Bifocal  Realy  Mirror  satellite  attitude  determination  plots 
%analyze  performance  of  dynamic  gyro  and  Kalman  filter 

%DG  momentum  error  or  gyro  bias 

if  dg==l 

figure(5) 

elf 

Hm=sqrt(H(:,l).A2+H(:,2).A2+H(:,3).A2); 

Hmdg=sqrt(Hdg(: ,  1).  A2+Hdg(:  ,2) .  A2+Hdg(:  ,3)  •  A2) ; 

plot(tout,Hm-Hmdg) 

title(' Angular  Momentum  Error') 

xlabel('Time  (sec)') 

ylabel('Nms') 

axis([0, stoptime, -.4, .4]) 

grid  on 

else 

figure(5) 

elf 

plot(tout,wb) 
title('Gyro  Biases') 
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legend(' wb  1' , '  wb2' , '  wb3 ') 

xlabel('Time  (sec)') 

ylabel('rad/s') 

grid  on 

end 

%KF  states 

figure(6) 

elf 

subplot(2 1  l),plot(tout,x(:,  1 :3)) 
title('State  Attitude  Errors') 
Iegend('xar,'xa2','xa3') 
xlabel('Time  (sec)') 
ylabel('rad') 

axis([0, stoptime, -.0008, .0008]) 
grid  on 

subplot(212),plot(tout,x(:,4:6)) 
title('State  Rate  "Bias"  Errors') 
Iegend('xbl','xb2','xb3') 
xlabel('Time  (sec)') 
ylabel('rad/s') 

axis([0, stoptime, -.00008, .00008]) 
grid  on 

%AD  errors 

figure(7) 

elf 

plot(tout,(q-qkf)) 
title('Quatemion  Errors') 
Iegend('qr,'q2','q3','q4') 
xlabel('Time  (sec)') 
axis([0, stoptime, -3e-3,3e-3]) 
grid  on 
figure(8) 
elf 

plot(tout,w-wkf) 
title('Rate  Errors') 
Iegend('wl','w2','w3') 
xlabel('Time  (sec)') 
ylabel('rad/s') 

axis([0, stoptime,- 1 .5e-3, 1 .5e-3]) 
grid  on 


3.  Attitude  Control  Performance  Results 

%Bifocal  Realy  Mirror  satellite  attitude  control  plots 
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%analyze  commanded  versus  controled  attitude 

figure(9) 

elf 

plot(tout,Mws) 
title('Wheel  Torques') 
xlabel('Time  (sec)') 
ylabel('Nm') 

legend('Wheel  1', 'Wheel  2','Wheel  3') 

grid  on 

figure(lO) 

elf 

plot(tout,hws) 

title('Wheel  Angular  Momentums') 
xlabel('Time  (sec)') 
ylabel('Nms') 

legend('Wheel  1', 'Wheel  2','Wheel  3') 

grid  on 

figure(ll) 

elf 

plot(tout,(q-qc)) 

title('Quaternion  Error') 

Iegend('ql','q2','q3','q4') 

xlabel('Time  (sec)') 

axis(  [0,  stoptime, -4e-3 ,4e-  3  ] ) 

grid  on 

figure(12) 

elf 

plot(tout,w-wc) 
title('Rate  Error') 
legend('w  1  ','w2','w3') 
xlabel('Time  (sec)') 
ylabel('rad/s') 

axis([0, stoptime,- 1 .5e-3, 1 .5e-3]) 
grid  on 
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